{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> <h1> Linear regression </h1></center>\n",
    "# 1. Model specification\n",
    "The linear regression model has the following form.\n",
    "$$\n",
    "p(y|\\vec{x},\\vec{\\theta})=\\mathcal{N}(y|\\vec{w}^T\\vec{x},\\sigma^2)\n",
    "$$\n",
    "Linear regression can be made to model non-linear relationships by replacing $\\vec{x}$ with some non-linear basis function $\\phi(\\vec{x})$, that is\n",
    "$$\n",
    "p(y|\\vec{x},\\vec{\\theta})=\\mathcal{N}(\\vec{w}^T\\phi(\\vec{x}),\\sigma^2)\n",
    "$$\n",
    "## 1.1 MLE estimation\n",
    "For training dataset $D=\\lbrace (\\vec{x}_i,y_i) \\rbrace_{i=1}^N$, the log likelihood function is \n",
    "\\begin{align*}\n",
    "log\\,p(D|\\vec{w},\\sigma^2) &=\\sum_{i=1}^N log \\left[ \\mathcal{N}(y_i|\\vec{w}^T\\vec{x}_i,\\sigma^2) \\right] \\\\\n",
    "                           &=\\sum_{i=1}^N log \\left[\\frac{1}{\\sqrt{2 \\pi \\sigma^2}} exp \\left[-\\frac{(y-\\vec{w}^T\\vec{x})^2}{2\\sigma^2} \\right] \\right] \\\\\n",
    "                           &=-\\frac{N}{2}log(2\\pi\\sigma^2)-\\frac{\\sum_{i=1}^N (y-\\vec{w}^T\\vec{x})^2}{2\\sigma^2}\n",
    "\\end{align*}\n",
    "We defined the **residual sum of squares** and **mean squared error (MSE)** as following\n",
    "\\begin{align*}\n",
    "RSS(\\vec{w}) &=\\sum_{i=1}^N(y-\\vec{w}^T\\vec{x})^2 \\\\\n",
    "MSE &=\\frac{RSS}{N}\n",
    "\\end{align*}\n",
    "The MLE estimation for $\\vec{w}$ is the one that minimizes the $RSS$, so this method is known as least squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHVCAYAAAA+QbhCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xl4VNXh//H3IXsCIRNIAAPIKpvIFhTFBVkU2dQqgvhVXFpsa22tWsWlgkv7q1pbq7XWlbqVTRFwR3BpVaRGQRZBNhXDlrBvCZDk/P44E0ggkEkyyczc+byeZ56ZuXPnzpmb5XPOueeea6y1iIiISHiqF+oCiIiIyLEpqEVERMKYglpERCSMKahFRETCmIJaREQkjCmoRUREwpiCWkREJIwpqEVERMKYglpERCSMxYa6AACNGze2rVq1CnUxRERE6sSXX365xVqbEci6YRHUrVq1IicnJ9TFEBERqRPGmB8CXVdd3yIiImFMQS0iIhLGFNQiIiJhTEEtIiISxhTUIiIiYSwsRn2LiISDXbt2kZeXx8GDB0NdFIlgcXFxZGZmkpqaGpTtKahFRHAhvXnzZrKyskhKSsIYE+oiSQSy1lJQUMD69esBghLW6voWEQHy8vLIysoiOTlZIS3VZowhOTmZrKws8vLygrJNBbWICHDw4EGSkpJCXQzxiKSkpKAdQlFQi4j4qSUtwRLM3yUFtYiISBhTUIuIiISxSoPaGNPBGLOozG2XMeYmY0y6MeZ9Y8wq/73Pv74xxjxmjFltjFlsjOlZ+19DRESmTZvGv/71r6Bu86OPPsIYw9KlS6v0vn/9618YY9izZ09Qy3Ms33//PcYY3nzzzTr5vLpUaVBba7+11na31nYHegH7gNeB8cA8a217YJ7/OcAFQHv/bRzwZG0UXEREyquNoO7Zsyfz58+nbdu2VXrf0KFDmT9/PsnJyUEtTzSq6nnUA4A11tofjDEXAv38y18APgJuBy4EXrTWWuBzY0yaMaaZtXZjkMosIiI1cPDgQerVq0dMTEyl66amptKnT58qf0ZGRgYZGQFdblkqUdVj1KOByf7HTUrD13+f6V+eBfxY5j25/mUiIlJLrr76al577TU+/vhjjDEYY5g4cSIA/fr149JLL+Xpp5+mbdu2JCYmsmHDBlasWMHo0aNp0aIFycnJdOnShUcffZSSkpJD262o69sYw9/+9jfuvPNOMjIyyMzM5IYbbmD//v2H1jmy67u0a3ratGlcf/31NGzYkObNmzNhwoRynwcwffp02rdvT1JSEueeey4LFy7EGFPl3oLi4mImTpxIy5YtSUhIoEuXLvz73/8ut86yZcsYPHgw6enppKSk0KlTJ5544olDr3/yySecddZZpKamkpqaSvfu3Zk+fXqVylFTAbeojTHxwAjgjspWrWCZrWB743Bd47Rs2TLQYoiISAV+//vfs27dOnbs2ME//vEPAJo3b37o9U8//ZQ1a9bw4IMPkpycTMOGDVm5ciUdOnTgiiuuoEGDBixatIgJEyZQUFDAHXcc/1/9I488Qv/+/Xn55ZdZvHgxd9xxByeeeCK33Xbbcd932223cckll/Dqq68yb9487rvvPrp06cJll10GQE5ODqNHj+bSSy/l8ccfZ/ny5YwaNapa++See+7hoYceYsKECfTu3ZvXXnuNK664AmMMl19+OQAjRoygY8eOvPzyyyQkJPDtt9+ya9cuwM1WN2zYMC688ELuuecerLUsWbKEHTt2VKs81VWVru8LgK+stZv9zzeXdmkbY5oBpVOw5AItyryvObDhyI1Za58GngbIzs4+KsirbcEC6NPH3Z96atA2KyJR6KabYNGi0Hx29+7w6KMBr962bVvS09MpKSmpsKt6x44dLFy4kKZNmx5aNmDAAAYMGAC4qS/PPPNM9u3bxzPPPFNpULdq1epQC/f888/n008/ZcaMGZUG9dlnn80jjzwCwKBBg3j33XeZMWPGoaB+8MEH6dSpE1OmTMEYw+DBgzl48CC33357wPsCYNu2bTz66KPcfffd3H333YfKmZuby8SJE7n88svZsmULa9euZebMmXTt2vXQPim1cuVKdu7cyd///ncaNGgAwHnnnVelcgRDVbq+L+dwtzfAbGCs//FYYFaZ5Vf5R3/3AXbW6fHpt9929++8U2cfKSIS7nr16lUupAEKCwuZMGEC7dq1IyEhgbi4OO666y6+++47ioqKjru9IwOrc+fO5ObmVlqOyt73xRdfMHz48HIThowYMaLS7R5p6dKl7Nu3j5EjR5ZbPmrUKFauXEleXh7p6em0aNGCn//850ydOvWoKT/btm1L/fr1GTNmDLNmzarzlnSpgFrUxphkYBBwfZnFfwKmGWOuA9YBpXvjbWAIsBo3QvyaoJVWRKQuVaFFG+6aNGly1LLbb7+dZ599lgkTJtCzZ0/S0tKYNWsWDzzwAIWFhdSvX/+Y20tLSyv3PD4+nsLCwkrLUdn7Nm3adNQgtOoMStu40bUPj/zepc+3b99OZmYmc+bM4a677uLaa6+loKCAvn378thjj9GjRw98Ph9z5szh3nvv5bLLLqOkpITzzjuPxx9/nDZt2lS5TNUVUIvaWrvPWtvIWruzzLKt1toB1tr2/vtt/uXWWnuDtbattbartTantgovIiKBqWhKy+nTp3PjjTdy2223MXDgQLKzs4mNDe1FFZs2bUp+fn65ZUc+D0SzZs0Ajmolb97sjt6mp6cD0LFjR1577TV27NjB3LlzKSwsZOjQoYcGuJ1++um8++677NixgxkzZrBy5UrGjBlT5fLUhGYmExHxiEBbtaUKCgpISEg49Ly4uJgpU6bURtEC1rt3b9544w3cGb7O7Nmzq7ydk08+meTk5KNGaE+bNo2TTjrpqFZ6XFwc/fv35+abb2bjxo1HdXMnJSUxfPhwrr32Wr755psql6cmdD1qERGP6NixI7NmzWLmzJk0b96cE044gRNOOOGY6w8aNIgnnniCdu3akZ6ezhNPPFHuFKtQuP322znttNMYPXo011xzDcuXL+eZZ54BoF69wNuW6enp3HTTTTzwwAPExsaSnZ3NjBkzePvtt5k82Q23Wrx4MbfeeiujRo2iTZs2bN++nQcffJBu3bqRnp7OW2+9xfPPP89FF11Ey5YtWb9+PU899RT9+/evle9+LApqERGP+OUvf8nChQu59tpr2b59OxMmTDh0LnVFHn/8cX7+859zww03kJSUxNixY7n44osZN25c3RX6CNnZ2UyePJk777yTWbNmkZ2dzZNPPsmgQYNITU2t0rbuu+8+YmNjefLJJ9m8eTPt2rXj5ZdfZvTo0YDrZm/SpAl/+MMf2LBhA2lpaZx77rk8+OCDALRr1w5jDHfeeSd5eXlkZGQwbNgw/vjHPwb9ex+PKdu9ECrZ2dk2JydIh7InTID77oOJE91jEZEALF++nE6dOoW6GFKBl19+mSuvvJK1a9fSunXrUBcnYMf7nTLGfGmtzQ5kO2pRi4hIWPnFL37BoEGD8Pl8fPXVVzzwwAMMHTo0okI6mBTUIiISVrZu3covf/lLtm7dSqNGjRg1ahQPPfRQqIsVMt4L6jDoyhcRkeqbNm1aqIsQVnR6loiISBjzXlBXcFK/iIhIpPJeUIuIiHiIglpERCSMeS+oNZhMREQ8xHtBLSIi4iEKahERkTCmoBYRkXImTpxI48aNDz3/6KOPMMawdOnS477v1ltvpVWrVlX6rLy8PCZOnMj3339fjZIeLdCyRhLvBbVOzxIRCaqePXsyf/582rZtG/Rt5+Xlce+99wYtqL3IezOTiYhIUKWmptKnT59QFyNqea9FLSIShSZNmkRCQgI7duwot3zZsmUYY5g3bx4Ab731FoMGDSIzM/NQAM+ZM+e4266oO3nHjh2MGTOGlJQUmjVrxh/+8Iej3rdx40auvfZa2rRpQ1JSEieddBJ33303Bw4cAOD777+na9euAJx77rkYYzBlekW3bdvG9ddfT5MmTUhMTOSMM85gwYIFVd43+/bt49e//jVNmzYlMTGR3r17H/WdP/nkE8466yxSU1NJTU2le/fuTJ8+/dDrs2fPplevXqSkpODz+TjttNP4+OOPq1yW6vBeUOv0LBGJQj/5yU8AeP3118stnzp1KpmZmfTr1w+A7777juHDh/PSSy/x2muvccYZZ3DBBRfw6aefVunzrrnmGt555x0effRRnn76aebMmcOUKVPKrbNlyxbS09P5y1/+wrvvvsvvfvc7Jk2axI033ghAs2bNeOWVVwB44oknmD9/PvPnzwdg//79DBw4kPfff5+HH36YmTNnkpGRwcCBA9m0aVOVyvqzn/2MSZMmcdddd/H666/TokULhg4dyieffALArl27GDZsGG3atOG1117j1Vdf5corrzxU6VmzZg2XXnop/fv354033uCVV15h2LBhbNu2rUrlqC51fYuIHMNN797Eok2LQvLZ3Zt259HBjwa8fsOGDRk8eDBTp07lmmuuObR86tSpjBw5kpiYGAB+9atfHXqtpKSEc889l2XLlvHcc8/Rt2/fgD5r2bJlzJw5kylTpjBq1CjAtYhbtmxJamrqofW6du3Kn//850PP+/btS0pKCtdeey2PP/44CQkJnHLKKQB07ty5XPf6yy+/zNKlS1m2bBnt27cHYODAgXTo0IFHHnmEhx9+OKCyLl++nMmTJzNp0iTGjh0LwPnnn88pp5zC/fffz3vvvcfKlSvZuXMnf//732nQoAEA55133qFtLFy4kAYNGpT7zCFDhgT0+cHgvRa1iEiUGjVqFPPmzWPLli0ALFq0iJUrVx4KU4Dc3FzGjh1LVlYWsbGxxMXFMWfOHFauXBnw53zxxRcAjBgx4tCy+vXrM2jQoHLrWWt59NFH6dy5M0lJScTFxXHFFVewf/9+1q1bd9zPmDt3Lr169aJ169YUFRVRVFQEwDnnnENOTk6VymqtZeTIkYeW1atXj5EjRx5qUbdt25b69eszZswYZs2addThg65du7Jz507Gjh3LnDlz2Lt3b8CfHwxqUYuIHENVWrThYMSIEcTFxTFjxgzGjRvH1KlTycrK4swzzwRcC3rEiBHs3r2b++67j3bt2pGSksI999xDXl5ewJ+zadMmGjRoQFJSUrnlmZmZ5Z4/+uij3HrrrYwfP55zzjkHn8/HF198wQ033EBhYeFxP2PLli18/vnnxMXFHfVaVUafb9y4kfr165OcnFxueZMmTdi3bx/79+/H5/MxZ84c7r33Xi677DJKSko477zzePzxx2nTpg0dOnRg1qxZ/OlPf2LIkCHExcVx8cUX87e//Y2MjIyAy1Jd3gtqnZ4lIlGqfv36DB06lKlTpzJu3DimTZvGZZdddmiA1urVq1m4cCHvvPMOgwcPPvS+goKCKn1O06ZN2b17NwUFBeXC+siwnz59OiNHjiw30Oybb74J6DPS09PJzs7mySefPOq1hISEgMvarFkz9uzZw759+8qF9ebNm0lOTj60rdNPP513332XgoIC5s6dy80338yYMWP4/PPPARg6dChDhw5l586dvPXWW9x0003ceOONRx2Xrw3e6/rWYDIRiWKjR4/m448/5o033mDt2rWMHj360GulgVw26H744YcqDyTr3bs34EZCl9qzZw/vv/9+ufUKCgqOCtXSwWOl4uPjAY5qYQ8YMIDVq1fTsmVLsrOzy91KR4oHWlZjDK+++uqhZdZaXn311UM9DWUlJSUxfPhwrr322gorFQ0bNmTMmDFcfPHFAVc6asp7LWoRkSg2dOhQkpOTuf7662ndujWnnnrqodc6duxI8+bNueWWW7j//vvZvXs3EyZMICsrq0qf0aVLF0aMGMEvfvELdu3aRbNmzXj44YeP6l4eNGgQjz32GKeddhpt27bllVdeYfXq1eXWadmyJUlJSbzwwgs0bNiQuLg4srOzueqqq/jnP/9Jv379uPXWW2nTpg1bt27lf//7H02bNuW3v/1tQGXt1KkTl19+Ob/61a/YtWsX7dq145lnnmHFihWHWutvvfUWzz//PBdddBEtW7Zk/fr1PPXUU/Tv3x+Ap556ivnz5zN48GBOOOEEVq1axfTp07nqqquqtN+qzVob8luvXr1s0Pz+99aCtRMnBm+bIuJ533zzTaiLEDRXXHGFBez48eOPeu1///uf7d27t01MTLTt2rWzkyZNsmPHjrVl/w9PmDDBNmrU6NDzDz/80AJ2yZIlh5Zt27bNjho1yiYnJ9vMzEx777332ltuucWeeOKJh9bZvXu3vfrqq63P57M+n89ed9119o033jhqWy+//LJt3769jYuLsy6WnB07dthf//rXtnnz5jYuLs5mZWXZiy++2H7yySfH/O4VlXXv3r32V7/6lc3MzLTx8fG2V69e9t133z30+ooVK+wll1ximzdvbuPj421WVpa9/vrr7datW6211n722Wd2yJAhtlmzZjYhIcG2atXK3nbbbbawsPC4P4fj/U4BOTbAjDQ2DLqKs7OzbVVG8R3XPffA/ffDxIkwYUJwtikinrd8+XI6deoU6mKIhxzvd8oY86W1NjuQ7XjvGLWIiIiHeDeoNfpbREQ8wLtBHQZd+iIiIjXl3aAWERHxAAW1iIhfOAyuFW8I5u+SglpEBIiLi6vyDF0ix1JQUFDh9KfV4d2g1mAyEamCzMxM1q9fz759+9Sylmqz1rJv3z7Wr19/1Nzn1eXdmcn0hyYiVVB6ecYNGzZw8ODBEJdGIllcXBxNmjQpd8nPmvBuUIuIVFFqamrQ/rmKBIt3u75FREQ8QEEtIiISxhTUIiIiYUxBLSIiEsa8G9Q6PUtERDzAu0Gt07NERMQDvBvUIiIiHqCgFhERCWMKahERkTDm3aDWYDIREfEA7wa1iIiIByioRUREwph3g1qnZ4mIiAcEFNTGmDRjzKvGmBXGmOXGmNONMenGmPeNMav89z7/usYY85gxZrUxZrExpmftfgURERHvCrRF/TfgXWttR6AbsBwYD8yz1rYH5vmfA1wAtPffxgFPBrXEgdJgMhER8YBKg9oYkwqcDTwHYK09YK3dAVwIvOBf7QXgIv/jC4EXrfM5kGaMaRb0kouIiESBQFrUbYB8YJIxZqEx5lljTArQxFq7EcB/n+lfPwv4scz7c/3LyjHGjDPG5BhjcvLz82v0JURERLwqkKCOBXoCT1prewB7OdzNXZGK+pyPGtllrX3aWpttrc3OyMgIqLAiIiLRJpCgzgVyrbUL/M9fxQX35tIubf99Xpn1W5R5f3NgQ3CKKyIiEl0qDWpr7SbgR2NMB/+iAcA3wGxgrH/ZWGCW//Fs4Cr/6O8+wM7SLvI6pdOzRETEA2IDXO9G4BVjTDywFrgGF/LTjDHXAeuAkf513waGAKuBff51687VV8P998N339Xpx4qIiNQGY8Og5ZmdnW1zcnKCs7F16+DEE93jMPhuIiIiRzLGfGmtzQ5kXe/NTKbzp0VExEO8F9QiIiIe4r2gVotaREQ8xHtBLSIi4iHeC2q1qEVExEMU1CIiImHMe0EtIiLiId4LarWoRUTEQ7wX1CIiIh7ivaBWi1pERDzEe0EtIiLiId4LarWoRUTEQxTUIiIiYcx7QS0iIuIh3gtqtahFRMRDvBfUIiIiHuK9oFaLWkREPMR7QS0iIuIh3gtqtahFRMRDFNQiIiJhzHtBLSIi4iHeC2q1qEVExEO8F9QiIiIe4r2gVotaREQ8xHtBLSIi4iHeC2q1qEVExEMU1CIiImHMe0EtIiLiId4LarWoRUTEQ7wX1CIiIh7ivaBWi1pERDxEQS0iIhLGvBfUZVkb6hKIiIjUiPeCWi1qERHxEO8FdVlqUYuISITzXlCXbVErqEVEJMJ5L6jLUlCLiEiE815Qq0UtIiIe4u2gFhERiXDeC+qy1KIWEZEI572gVte3iIh4iPeCuiwFtYiIRDjvBbWOUYuIiId4L6jLUotaREQinPeCWseoRUTEQxTUIiIiYcx7QV2WglpERCKc94Jag8lERMRDAgpqY8z3xpglxphFxpgc/7J0Y8z7xphV/nuff7kxxjxmjFltjFlsjOlZm1/guNSiFhGRCFeVFvW51tru1tps//PxwDxrbXtgnv85wAVAe/9tHPBksAobEB2jFhERD6lJ1/eFwAv+xy8AF5VZ/qJ1PgfSjDHNavA51aegFhGRCBdoUFtgjjHmS2PMOP+yJtbajQD++0z/8izgxzLvzfUvK8cYM84Yk2OMycnPz69e6SuiY9QiIuIhsQGu19dau8EYkwm8b4xZcZx1K0rKo5q21tqngacBsrOza6fpqxa1iIhEuIBa1NbaDf77POB14FRgc2mXtv8+z796LtCizNubAxuCVeAqUVCLiEiEqzSojTEpxpgGpY+B84ClwGxgrH+1scAs/+PZwFX+0d99gJ2lXeR1TkEtIiIRLpCu7ybA68Yd+40F/m2tfdcY8wUwzRhzHbAOGOlf/21gCLAa2AdcE/RSi4iIRIlKg9pauxboVsHyrcCACpZb4IaglK6m1KIWEZEI572ZycpSUIuISITzdlD/+GPl64iIiIQxbwd1z9DNXioiIhIM3g5qERGRCKegFhERCWMKahERkTDmzaAeMSLUJRAREQkKbwb1kCGhLoGIiEhQeDOoRUREPMKbQa1LXYqIiEcoqEVERMKYglpERCSMKahFRETCmDeDWkRExCO8GdRqUYuIiEcoqEVERMKYglpERCSMeTOoRUREPMKbQa0WtYiIeISCWkREJIwpqEVERMKYN4NaRETEIzwX1JOXTCbrh9+wOSXUJREREak5zwU1wIbi7WxLCnUpREREas5zQe1L8gGwIzHEBREREQkC7wV1ogvq7WpRi4iIB3guqNMS0wC1qEVExBs8F9SlXd/bFdQiIuIBngvq0ha1ur5FRMQLPBfU8THxJJfEqOtbREQ8wXNBDeArjlfXt4iIeIIngzqtJF5d3yIi4gmeDGqfTVDXt4iIeIJHgzpRXd8iIuIJngzqNJtwuOvb2pCWRUREpCY8GdS+kvjDXd8KahERiWCeDeqdiVBsUFCLiEhE82RQpxXHA7ArIcQFERERqSFPBrWvKBbwz06mFrWIiEQwTwZ1WnEcoPm+RUQk8nkyqH3d+gD+K2ipRS0iIhHMm0E99FJAF+YQEZHI58mgTktsCPi7vtWiFhGRCObJoPYlumtSq+tbREQinSeDOiUuhdhif9f3kiWhLo6IiEi1eTKoTb16pBX6W9S9e4e6OCIiItXmyaAG8BXq9CwREYl8ng3qtEKN+hYRkcgXcFAbY2KMMQuNMW/6n7c2xiwwxqwyxkw1xsT7lyf4n6/2v96qdop+fL4CdE1qERGJeFVpUf8GWF7m+YPAX6217YHtwHX+5dcB26217YC/+terc+r6FhERLwgoqI0xzYGhwLP+5wboD7zqX+UF4CL/4wv9z/G/PsC/ft0xRl3fIiLiCYG2qB8FbgNK/M8bATustUX+57lAlv9xFvAjgP/1nf7161Rp17fOohYRkUhWaVAbY4YBedbaL8surmBVG8BrZbc7zhiTY4zJyc/PD6iwVeErhIMxsC8u6JsWERGpM4G0qPsCI4wx3wNTcF3ejwJpxphY/zrNgQ3+x7lACwD/6w2BbUdu1Fr7tLU221qbnZGRUaMvUZG0QnevAWUiIhLJKg1qa+0d1trm1tpWwGjgA2vtFcCHwKX+1cYCs/yPZ/uf43/9A2vreB5PY/AVuIc6Ti0iIpGsJudR3w7cbIxZjTsG/Zx/+XNAI//ym4HxNSti9ZS2qDXyW0REIlls5ascZq39CPjI/3gtcGoF6xQCI4NQthrxqetbREQ8wLMzk6nrW0REvMCzQa2ubxER8QJvBrV/whNQ17eIiEQ2bwY1EGMhVbOTiYhIhPNsUIP/ClpqUYuISATzdFD7CtX1LSIikc3TQa0Lc4iISKTzZlD7L9ala1KLiEik82ZQ++ma1CIiEuk8HdTq+hYRkUjn6aD2FcDeeDhYfDDURREREakWbwd16aQnhTtCWxAREZFq8mZQ+weTHZpGtHB7CAsjIiJSfd4Mar/SC3OoRS0iIpHK00F9qEVdoBa1iIhEJk8HtY5Ri4hIpPN2UJdek1rHqEVEJEJ5M6iPHEymrm8REYlQ3gxqv6QiSChS17eIiEQuTwc1uO5vdX2LiEik8nxQpxXC9j35oS6GiIhItXg+qH2FsGPBx6EuhoiISLXEhroAtS2tEPJi1PUtIiKRyfst6gJd6lJERCKX94O6EHYoqEVEJEJ5PqjT/EFdYktCXRQREZEq83xQ+wqgpB7sObAn1EURERGpMu8HtWYnExGRCOb5oNY1qUVEJJJ5Pqh1TWoREYlkng9qXZhDREQimeeD2qeubxERiWDeD2p1fYuISATzfFA3OADGqutbREQik+eDup71T3qiFrWIiEQgzwc16JrUIiISuaIiqNMKFdQiIhKZoiKofer6FhGRCBUVQZ1WqMFkIiISmaIiqHWMWkREIlV0BLW6vkVEJEJFRVCnFUJhUSGFRYWhLoqIiEiVREVQa3YyERGJVNER1Lowh4iIRKioCGpdk1pERCJVVAS1ur5FRCRSRUVQ65rUIiISqaIiqHVNahERiVRREdSlLWp1fYuISKSpNKiNMYnGmP8ZY742xiwzxtzrX97aGLPAGLPKGDPVGBPvX57gf77a/3qr2v0KlYsvhmSToK5vERGJOIG0qPcD/a213YDuwGBjTB/gQeCv1tr2wHbgOv/61wHbrbXtgL/61ws5X0yyur5FRCTiVBrU1tnjfxrnv1mgP/Cqf/kLwEX+xxf6n+N/fYAxxgStxNXkq5eirm8REYk4AR2jNsbEGGMWAXnA+8AaYIe1tsi/Si6Q5X+cBfwI4H99J9Cogm2OM8bkGGNy8vPza/YtApBWL0UtahERiTgBBbW1ttha2x1oDpwKdKpoNf99Ra1ne9QCa5+21mZba7MzMjICLW+1+eolq0UtIiIRp0qjvq21O4CPgD5AmjEm1v9Sc2CD/3Eu0ALA/3pDYFswClsTaTEpGkwmIiIRJ5BR3xnGmDT/4yRgILAc+BC41L/mMNKpAAAgAElEQVTaWGCW//Fs/3P8r39grT2qRV3XfPU0mExERCJPbOWr0Ax4wRgTgwv2adbaN40x3wBTjDEPAAuB5/zrPwe8ZIxZjWtJj66FcleZr14Ku/bvorikmJh6MaEujoiISEAqDWpr7WKgRwXL1+KOVx+5vBAYGZTSBVFavWQAdu7fSXpSeohLIyIiEpiomJkMXIsaNN+3iIhElugJ6hgX1Br5LSIikSRqgrq061sDykREJJJETVD7/EGtFrWIiESSqAnqQy1qHaMWEZEIEjVB7TPq+hYRkcgTNUGdYhKIrRerrm8REYkoURPUxhjSEtPU9S0iIhElaoIawJfoU9e3iIhElOgJamvxJfnU9S0iIhEleoJ60ybX9a0WtYiIRJDoCerx4/EVx6tFLSIiESV6ghpIO1BPg8lERCSiRFVQ+2Lrs71wO2FweWwREZGARFdQx9SnqKSIfQf3hbooIiIiAYmqoE7zX0FLA8pERCRSRFVQ+2LrA5rvW0REIkd0BXWMC2qN/BYRkUgRVUGtrm8REYk0URXUvtKgVte3iIhEiKgK6rR6LqjV9S0iIpEiuoI6RtekFhGRyBJVQR1DPVITUtWiFhGRiBFVQU1JiS7MISIiESW6gtpad01qDSYTEZEIEV1BvWmTrkktIiIRJbqC+mc/U9e3iIhElOgKalDXt4iIRJSoC+q0xDR1fYuISMSIuqD2JfrYe3AvB4sPhrooIiIilYq+oE7yAZqdTEREIkPUBXVaYhqg2clERCQyRF1Q+xJdi1oDykREJBJEX1Cr61tERCJI1AW1ur5FRCSSRF1Qq+tbREQiSdQFdWmLWl3fIiISCaIuqJPikkiISVDXt4iIRISoC2pAF+YQEZGIEZVBrQtziIhIpIjKoNaFOUREJFJEX1B//bW6vkVEJGJEX1B3766ubxERiRjRF9So61tERCJHVAZ1WmIaO/fvpMSWhLooIiIixxWVQe1L9FFiS9i9f3eoiyIiInJc0RnU/gtz6Di1iIiEu6gMak0jKiIikaLSoDbGtDDGfGiMWW6MWWaM+Y1/ebox5n1jzCr/vc+/3BhjHjPGrDbGLDbG9KztL1FVujCHiIhEikBa1EXALdbaTkAf4AZjTGdgPDDPWtsemOd/DnAB0N5/Gwc8GfRS15DvzbmAWtQiIhL+Kg1qa+1Ga+1X/se7geVAFnAh8IJ/tReAi/yPLwRetM7nQJoxplnQS14DaS9MAXSMWkREwl+VjlEbY1oBPYAFQBNr7UZwYQ5k+lfLAn4s87Zc/7IjtzXOGJNjjMnJz8+veslrwHcgBlDXt4iIhL+Ag9oYUx94DbjJWrvreKtWsMwetcDap6212dba7IyMjECLERQNimIwGHV9i4hI2AsoqI0xcbiQfsVaO8O/eHNpl7b/Ps+/PBdoUebtzYENwSlucNTDaBpRERGJCIGM+jbAc8Bya+1fyrw0GxjrfzwWmFVm+VX+0d99gJ2lXeThxJfkU1CLiEjYiw1gnb7AlcASY8wi/7I7gT8B04wx1wHrgJH+194GhgCrgX3ANUEtcZCkJaap61tERMJepUFtrf2Eio87AwyoYH0L3FDDctUuY3RhDhERiQhROTMZoGtSi4hIRIjOoC4pIS1Bg8lERCT8RWdQr1jhBpOp61tERMJcdAY1br7v/cX7KSwqDHVRREREjilqgzqtSLOTiYhI+IvaoPZN0nzfIiIS/qI2qNNWuenINfJbRETCWdQGte/HLYC6vkVEJLxFb1D7x5Cp61tERMJZ1AZ1mj+o1fUtIiLhLGqD2lfg7tX1LSIi4SxqgzquBFJik9WiFhGRsBa1QQ2QltBQx6hFRCSsRXVQ+zTft4iIhLmoDuq0hFR1fYuISFiL6qD2JaRpMJmIiIS16A7qJavZnr8u1MUQERE5pqgO6rRla9ixTy1qEREJX1Ed1L5C2JUIxSXFoS6KiIhIhaI7qP2TnuzcvzO0BRERETmGqA7q0mlENaBMRETCVVQHtS7MISIi4S6qg/rQhTkeeQB26HxqEREJP1Ed1IcuzPHuLJgwIbSFERERqUB0B3Vp13cSUKyR3yIiEn6iOqgPdX0nhrYcIiIixxLVQZ1yAGKLYXsisGABPPlkqIskIiJSTlQHtcF1f+9IBHJy4Je/DHWRREREyonqoAbX/b09KdSlEBERqVjUB7WvwN/1XSo+HiZODFVxREREyon6oE4rPGIw2cGDcO+9sG1byMokIiJSyrtBPWMGDBtW6Wq+Y3V9N2oEy5YFv1wiIiJV4N2gvvhiOPnkSlc7quu7rLPPdtsREREJEe8GdYBKu75tRS9u2wYzZ4Ix8NBD8NFHdVw6ERGJdrGhLkCo+QqhKAb2xkP9A8dZ8fbb3X3HjtC3L7RoAddcA3v3QuvWEBcHMTF1UmYREYkeCurS+b4TKwnqUitWuBuUHx3epAls3gy/+AUkJcGBAzB8uDvO3akTZGbCN9+47vT8fDe6vHlz2LABsrLAWigpgbQ0N6AtUdOliYiIgrrcNKItdtVgQ5s3u/uys5v9/e812CCQnu6637t3d4F+4AAMGgTvvw8jR7oW/PLlcMkl8PXX7ph88+aui374cNi0CQoLoUcPN/Na795uuz/8AKedBt995z6jcWNYuxY6d4Zd/p3QpAls3epeB1eRSElxlQn1HIiI1JmoD+pyF+YIN6WniC1adHjZ9Onu/plnDi/7+OOj3/vEE7VXLnAVgtxcF/5btrjbsGEwdy6MGOF6FRYtgssuc70KrVtD27bw4YcwdKhbf+dO9/6cHDjlFEhIgJUr3aGF776D1FRo1swt69LFHWY4eND1QOTnu0pETIyrwDRs6F6Lj3djCoyp3e8vIlJHoj6odWGOasrNdfdffHF42eTJ7v655w4v++STo9/7j3/UXrkAfD7Yvh3at4c9e2DjRhgyxFUihgxxhxfmz4cxY+Dbb6FpU+jaFf7zH/f6rl3uPWee6b5fly7QoIGrePTv73okkpLgxBNhyRLo1g0KClxFok0b995GjVylYfdut/19+9zz+HjXOxEXV7v7QEQ8w9tBbSscy11O2WPU4hHbt7v7VasOL3v7bXc/c+bhZRVdg3zSpNorV1lJSbB/vzuUcMopsHixO0SRkQHz5sHVV7vDESkp7jDF3LnuEMeBA7BuHQwY4HoiOnRwhy7mz4cLLnCViHr13KDHnBzIznbvyc93lZHvv3ef0aABrF/vKjM7d0JsrKvg7NjhKjLWul6JhAT1ToiEmLeDOgBh3fUt3lVQcPjx4sXufuHCw8vK9kq8+aa7L3uI47HHjt7mvfcGr3ylEhPdOIfYWHfI4Ycf4PTTXY/K/v2uh+G99+AnP3EVhGXLXIXiyy9dxaBFC1fu4cNdxaCwEHr1gk8/hT59XIVgzRp3uGPlSleJyMyEpUvdYZGtW11lpmVLV0E54QRXcdi71x0W2b3blTE+HoqKXAXIWlcWEY+I+qBuqK5vkWMr9P+BFBW5kAbXei81ZYq7L1ux+Oyzo7dTdkxFbYqJgeJi97h1azfWoXNn1zOwcCFceqmrRHToACedBHPmuEMg27a5Qx79+rnyn366e8+XX7qeiuXLXY9D69auktG/v6tE7NnjKiTLlrkxGHFxrhJzyinuEEhyshtLsXGjG9exf7+rSKSlucpGcrKrVBjjKkMiFTA2gO7h2padnW1zcnKCv+Hx4+HBBytdreF4uHoR/O3d4BdBRCQgpWMrwJ3SuXy562nYs8cdBrnoIvjgA3fYIy3NVSguv9yNs2jUyJ318eGHMHiw62nIzYWzznLjLDp3doc7vvrq8DiLhARX8Vi4EHr2dOModu1yFZh169whlaQkd9ikVSv3WmysG+S5e7e7L+29iI3VIZIqMsZ8aa3NDmRdVeE4znzfIiJ1pTSkwYU0uNZ7qZdecvcvvnh4WUUNnGefDX7ZKpOS4noIkpLcoYsffnCHNtascSHety+89ZY7C6S42M1FcdFF8Pnn7lBI06au4nHhhW4cRUmJOy31o4/gnHPcOIu1a90Az2XL3GGYxo3dAM++fd3psSUlbjDnihXQrp17vnOnq4zk5bmKSlKSq3BkZLjeopgYt+zgQVdxgbCscHi7RX377W7qz0p0/zm03AmzJwe/CCIiEiGMcb0E7drB6tVuTMSJJ7rDPYsXu8McQfuowFvUGnFBJRfmEBGR6FDacF292t1v2HB4TMa0aaEpEwpqoIJrUouIiJQVwjMJFNToGLWIiFQihMeuKw1qY8zzxpg8Y8zSMsvSjTHvG2NW+e99/uXGGPOYMWa1MWaxMaZnbRY+WNT1LSIixxXOQQ38Cxh8xLLxwDxrbXtgnv85wAVAe/9tHPAkESCtEPbFwwFda0JERCoSzl3f1tr/ANuOWHwh8IL/8QvARWWWv2idz4E0Y0yzYBW2tvg06YmIiBxPmLeoK9LEWrsRwH+f6V+eBfxYZr1c/7LQCPDUM833LSIi4SrYbfmKqhwVpqUxZpwxJscYk5Ofnx/kYlSNrqAlIiLHFc5d38ewubRL23+f51+eC7Qos15zYENFG7DWPm2tzbbWZmdkZFSzGMGhC3OIiMhxRWDX92xgrP/xWGBWmeVX+Ud/9wF2lnaRhzO1qEVE5LhCGNSVzvVtjJkM9AMaG2NygQnAn4BpxpjrgHXASP/qbwNDgNXAPuCaWihz0EXCMWoLbEmGjH2hLomISBQqKQnZR1ca1Nbay4/x0oAK1rXADTUtVNAEOpgsjLu+f0yFl0+BF7vBigyYNNNd6SvaWaC4HsSG7m9HRKLJwYMh+2hdPQtILIKEovDp+t4dDzM6uXD+sDVYA2f+AKflws+HQed8OHV9qEtZ9yywsBlMPhmmnuxCesEz6mXYFwfvt4HsDZC1O9SlEfGoAwdC9tEKaj9fAWysH7rPLzYwr40L5xmdoCAO2m6DiR/B/y2GNtthaxJkj4OfjIKcp6HpntCVty4tbwxTTna3lY0hthgGrnWVmMtGwvsvRV/LuqgezGsNr5wCr3eEPQnQOQ8WPAv1Q/f/JOQO1oO5bWBBc7jlM2gQxftCgkwt6tA7PRde7ua6wR9637Wy68LSTBfOr3SFDamQVgBXfe1up/9Y/ny3RgUwcwqc/lMYORLmvQjxxXVTzrr2XZprNU8+GRY3BWPh3O/g1s/gJ8vdvnixG4y9GH43CP76XqhLXPss8L8sF85Tu0Befff7MnopdNsMv74Arh8GL8+o+DxJryox8ElL97syvQtsTXbLl2XAtOnRtS+O9E0GvNYJkorc345Abqr7f1Ll3qf6oWvJKaj9Jr8K4wfCo6fDf06EqdOhw9ba+az8ZPfP9sVuris3thiGrIIr34VhK49fSei2GZ6fBZdfCjcNhn+8VTtlDIUNDWB6Z9dy/tx/kt/pP8Lf3oGRy6DZET0IV30NXzZzP7NeG13PgxetbOQqcv/uCqsbucM0w7+FMUvc702Cv7K2IxF+3x/6/gi//CK0Za5tFljU1O2TKSdDbkNIOggjvoXLl8CyTLhrAPztR7jp81CXtu5YYHETeLWzu60oc+ZregFcuzBkRQup3FS3P6Z3hs9aut6WP8+p4kY6d66VsgVCQe2XUOxaZQO+g6svgp7Xw9/fdgO3glUj3xcHj5wOD54Je+Mhez089rZrEVXlOOvopfBVM3i4L/TcCD/9KkgFDIH9MW6w3CunwEet3PH47hvhT+/DqGXQasfx3//nOe4f9s+GQ5c86LGpTopd6zbWdwH0766Qk+VaAP2/gzv/63oUGu4/+j13/hfmN3cVuOwN3hzHsLKRazlP7grf+g+DnL8G/jQXLvz2cLf/iG/hixNcb0vv9a7y4lUWyDkBXvOH85p0qFcC5/wAN/7P7YurL4Ibhrj/F9098jdSmfUN3P6Y1sWFM0C3TfDAPPe/pcoCHJxcG4wN4YeXys7Otjk5OcHf8O9+B3/+c5Xftr4B/N9P4KPWMGYxPPkWpFbwjzFQxQZe6gZ39Xfd25d84449n5xX6VuPu80hV7hw+/hf0Ce3+tsKhdJ9MrEf/JAGJ22By5e6SkjHLVXbVl4K9BoHMdYdu28cwYPL5raBB/vCB62hpB703ABXLHH75YQAuuq2JblKpgW+esodIoh0uamuq39yV/jyBFdpOfsH16NwyTfH/o47Et2YjoJYWPgUZO6t23LXphIDnzd3QfRaJ1iX5iot/b+DS7+Bi1aUr/znpbjfi8Qi9zdSOneE15SG8/Qu8Kk/nE/ZBJctg5HfwEk16SWdOhUuuywo5QQwxnxprc0OaF1PB/Wtt8Ijj1TrrcUG/t9ZMKEftN4BU151rZSqmtcabj0PFjWDU3PhkTlw5rpqFeko25Kg98/cwLMvnzq6azgcWWBmR9ctuTzD9Sr8cZ4bHFaTnosvToCzroW+6+C9lyNvcNl3aXDL+fB6J2ixE8YucgFd1UoLuNZV32td79Cb/4Z6of8Tr5ZV6XDz+fDWSa6npdcGF86jlgZ+fPHrJtDnp24MyhwPDDr8pKWrtMzo5Cr98UVw3hq4ZLlrOacfp2L2WQs452p3eG3GVO8cu1/fwPUmTCsTzl03+8N5WRAPYU6ZAqNGBWljVQtqrLUhv/Xq1cvWiltusdZ1WFT79t+W2Ba/xcb9Hvvn07HFJrD3Lc3ADhmDZSL2xJuwk0/GltSwLBXdFmdik+/Enn4dtjAm+NsP5m1ea+ypP3X7pOMN2Fc7BXefTOrutn3LeaH/roHe9sZhf38uNuFu93P845nYgtiab/fJbLcv7js79N+xOvvkrv7Y+LuxqeOx9/TDftuo5r8XdwwI/Xer7m1jfewll7nvkXgX9qJR2Je7YnckVG07f+3jtvHwGaH/TjW9vdcWe+Y1WDPBfaeuv3C/78sb19JnTp4c1HgCcqwNLCMDWqm2b+Ec1Bbs1iTsxaPcL8PgK7CbU4697qYU7PXDsPXuwTYc7/4ggvGP93i3aZ1d2cYNq93Pqe7tfydgB17pytjit9jnemAP1qudz7phiPucf58c+u99vFuJ/+fW4reuvGN+gv0xNbjb/7+L3T+xOW1C/30DLfPrHbEtb3L75MqLXUAFY9s/He62Ofuk0H/Pqu6T53pg0253lbk/nondHV+z7V06EhtzD/bjE0P//apz2xt3+O+87a9rOZzL3l55JajxpKAuFaSgLv0Ff6K3+2Npegt2buujf3keOAtb/w5s7O+xvx6MzU+uu1/eOwa4X9x/9qq7z6zs9k1j7E/8rYDGv3O1+dqutByohz3rGmzSXdhFTUK/Dyq6Lc7E9hvr9ku3n2P/07J2PmdPHLbLL92+XxfESkBt3FalYy+4wu2Tk38R/BApiMX2uN4F3hpf6L9vILc1PuyAq9w+OfvqmvUqlL3tTMC2v9H9HwtWRaiubjnNXG8cE7E3nV/7/0/K3TIzgxpPCupSQQzq0tvXTbCdbnAtlTsGuO7mF7phm/tbRhePwq5Mr/tf4CLj/tHF/R77SYu6//yytx8aYq+50PUqNLgDO/Ec7K4atAKqetuUgs26GdvqN9gtSaHdF2VvW5Owv7rA7Zf021z3dFGAh1Kqe1vRyFUe+1yH3R+Gh0ZKu/7j73a/K3/t4ypbtfFZa3wuqHtcX8f/4Kt4O1jPHWZLust1/f+zV+CH3AK9Lc502+83tvZ6t4K9Tx44yzWCsm4+uqFUZ7cgUlCXuvnmWvlh7YnDXjfCBXPD8e6+989qr2UU6G1bIradv6ac26B2PqPIuH+uW5Ow6xu4f37LMrBfNsN+2sLVcuPvdj0PN59Xt70KZW+fZ7lyDLwy9P+Iioz7Z9voNhfSNwxx+6+uPr/00MivB4d2P5S9lYCd1cFVppiIveIn2A110Lp74yT3eT8dHvp9UNFtURNs9s9cGUeMrr2/Y4trYETCsfvVPuwZ17qyjr7E/Z8LWXmCqCpBrfOoqyHlIDw7Gwatgb+f6iaXGLUs9KNrfYVu5rI+P4VLRsHHkw5PhhGI7YnudK95bdz5uLsSoDC2/K0o5vjbqFcC1yyCCR9Bi101+TY1c9p6NxnMTy90p8U9ODc05fi0Bdw4xE1sc8738Ng7cMrmui3DyG/gN5/D3/rAGT9W8xzSIFrjg99c4EZzd8mDjya5c37rwrCVcMd/3RkdZ/zoflfDQWEs3H82PNTXjdyeOt2NWK7NkdlXfe1Gkf+/s9zEQsNX1uKHVYMFnu/h5gWIsfDKa27UfzRSUNfAqGWh/6d3pC758MJMF9S/HOoqFMf6Y98X5/5Q57V24fxVM7AGkg+4f2Idt7jzLsveEoqPXlb21n6rO50tHFy30J13+9CZbuayy6r4syqMdae0zGvtTiUry5SplB25f0tf25EIc9tC850wZbr7/FCdEvPQ+2760Z+OcLPbVee0r5oqiHWT/fzpTIgrhkfegxsXQFwdnzJ134ewIMv9ffTYFPoJQP7bEn42wk3gcvVCN4lPXZ3//tg77nS+qy52592Hy99ufjKMGw4zO7mpg194PbQV/1Dz9nnUt9wCf/lL8LcbAX5/LjxwDjzx1uHpJA/Wc/+sP2h9uNV8INb90+yTCwPWunNvT13vnTnED8TAuWPd7GWfPwtdjzPJTLFx681t4/bPf1tCYRzElLiJEkrPwT3yL8aWSd+yrxng4uUw/hPXCxNquanQ43o38cf/nql6mVanw5snwTvtID/F9SCV3ow94jlHv7aisZuY4/IlLowCmcCltuSluH2RFIQJQLYmue+XcrBqfze7Ety0xU/2hlbb4ek3YNDa6pejutb63IRBbbbDp8/X3XUOjuWt9nDdha6i+//mwm8WhL638pAg5qUmPCkVxUFdYmDE5fBeW7j9U1jY1M1hvifB/ePsvulwMJ/1Q3gESW3ZWN/NUJVYBF88c3hSCIubbnFuG9dq/qA1bPNf0OHkzW7fDFzrZsGqycx04WRuGzjvSteF+FIlF+84WM/1uLx5kuum/raxW94xH9ptcxWUkjI3S/nnJab8OvUPuGlO+39XF9+0cp+2gH5Xw9BV8PqUwHs7io27OtfsDvDGSfBN5uHXYovd31LKAfd9Sx8fuSz5oJu4ZH2qm4v8/g9C+zc4uwNceDlcnwP/fLP62zkQA5vqu6sR1j9QtR6kvXFucqh/9nazib0yo2azN9YKBXUtBPX69W7Kt8+i87IxOxPg1J+5S0OetMUFz4C10O97b0wtWRWfN4ezr3EhMXbR4VbzD2nu9RY7XSgPWOvWiYRZ3qrrgbPdxTv+8Sb84og/u/xkeKe9C+f32sKuRDf7Vb/v3fHdoatcy8srHu0Dvx0MD74Pt3167PV2x8OctvBGB9fi25LiQvnsH9xc4/HFLmj2xh++3xN//GUdt8CTb7rxFOFg/EB3aOLFGXBlFS5wsy0J3m7vKi3vtnO/M+D2T3qBuzUqOPy4otuBGDcL3ep0d5Wv+z+o2viaOqOgroWgLmW8Mlle1e1McP8cqnxJNw96pieMG+EepxW4QB7o71Vov9U7UypWpsTAsDGuF+G/k9zVuEpbzZ83d63gprtdKA9b6faRV69xbYFRI9182fNedBWSUusauvB5owN82ModJvIVwAWr3HSd56/x1pzZRfVg0JWut2DBM8c/TLQq3d+j0MH1uhTXgyZ73FXdsje4bv1tSce+7Uk4epstd8CLr9fdwMJqUVArqKX2vd3eHaPtsdGNJI1WW/0X78hNdRf/ADfv+rCV7tZjUxgdF6xluxLcnPk7E11Q/LelC6Cvm7rX2291ATTiW3cVrkifL/x4NtV3x+4bHHDH7ksP+RQbmN/CVVxmdzh8+cyum91+GeEP6EB/Zw7EuLNMSoN7V4Lbt2F/iElBraAWqUuLmsJf+7gWzAWrvN3dX5mlmXDaT2FfvDvF8Mx17nSl4d/W3nXpw9V/ToT+Y91AyMuXumAu293f73sXzMNXVn4ZWk858UT4/vugbU5BfSQFtYhU4tMW8H0aDF4dfWM4jvTwGXDbee6xrwCGlHb3r674WuhRoXVrWBu8YflVCWqdRy0igut67ftjqEsRHm79DJrshRN3eL+7P2AhbNQqqEVEpByDm7lMyigJXW2lXsg+WUREJFIoqEVERMKYur49qG1b2LwZ7rkHMjPdIIQRI2DdOvD53AjCpUvhlFNg717YudO9Z+1aSE+HxERYtQo6d3YTtxQXQ1YWfPEFdOoEe/a4bZ18Mnz1lfuMBg3gv/+FM8+ElSuhsBDat4d33oE+fWDXLli8GM44Az7+GJo0gfr1YcYMGDgQliyBLVugSxd45RU45xxXrg8+gP793X1mJqSmwurV7nFeuE0dJCJSC0LYotao72Dp2BGeeQbi412QpqRAvSjpsCgqgpgY2LfPfWdrXUWiYUPYutUtS0qCH36A5s3dsoICaNbMVQ6aN3fbWLvWVRKWLXP7r3Fj+M9/4PTT3Xu3bXMVl7lzoVs3957PP4ezz3b36emQkeEqHv36ue2tWQM9e8K0aXDqqXDgALz5Jpx/Prz3nqusNGkC8+fDSSe5z9m/Hxo1cuUUEQH3f2JT8K7gotOzjlSbQb1smWtZNm5ce58hda/0CrSlteh9+yA52VVASkpcxWPTJlc52LvXLW/SxFUMfD5XOVm71lXg1q51z0srBD16QH4+bNzoKh6ffeZ6U2Jj3eOzznI9H7GxruflzTddxWPzZlixAnr1grffdtuOiYHZs13FY/58V+bWreGll2DwYFfGr76Cvn3h009dGWJjXS9N48auB0VEKpeZ6f4Gg0RBfaRgB/VZZ8Gzz7oWmEikKypyFYl9+1yIFxe7x6mpLsjj493yDRtc70denlsnI8MdvmnVylVWNm50fxNLlrgeifr1YcECOO00V4E5cADatHE9Ir17w+7d8M037vVPPoEWLdxnvvceDBoEy5e7Qy8nnQSvveYO2eza5Q7vnHOOq6y0bu3eM2UKDBvmPnvTJrAhoqQAAAkrSURBVNeL8tlnrmx79rjvkZHhKkgi1dG4cVB/fxTURwpWUA8c6EL6nnuCsz0RqXtFRe5/QmGhq4Ts3w8HD7rDLfn57r642B1qadrUVUDi4lzF49tv3biPvDxXOWnZEhYudBWCoiI3NqRbN9fTVr++++f+8ceukvHjj+UP3/To4d6zYIGreCxY4HpjmjRxvSTnngvffeduPXq4Qzq9e7v3vP2260X54AP3GU2buvecfbYbu/L99278y+LFkJbmeoF27XJl2hPFU9DVRHp6UA+HKaiPFKygzs9XF7eIRC5rXSWkpMQdNtmzxw1cPXDAVVxSU10lJDnZ/d/cvNkNYt20yT1PT3e9IG3auN6Obdtcr8aSJa6yEBfnHvfs6Sot9eq5XpgPP3Q9J1u2uApL9+5u/En79q6y9PHHbsDq4sXuPW3awBtvuIpHfj58/bUbEPvOO9Chg/ucmTNhyBA3wLaoyPW8vPKKq8Bs2QJffnn4kE/Tpu5w1XffwQknuN4hcN+9MIArqyQluXW3bQvaj0JBfaSaBnWjRpCb635QIiIS+UpKXDYcOOAqByUlLrSTk13vQ0yMqxDk57vlGze6Xo4g0RSiwXT77fB//6eQFhHxktKzchLKXHOz9HGjRoeXpaS4+44d66ZcFVBQV+ZPfwp1CUREJIopqI/lootg5MhQl0JERKJclMzIUQ39+sGYMaEuhYiIRDkF9bEMHx7qEoiIiERJUH/yCTz3XODrW+tODxAREQmx6Ajqvn3hsstCXQoREZEq02CysqZMCekQfBERkSMpqMsaNSrUJRARESknOrq+RUREIlT0BPXxpkrt0QPeeqvuyiIiIhKg6Anq42nRwk3uLiIiEmYU1ADt2oW6BCIiIhWKnsFkscf4qu+84y6vJiIiEoaip0WdlATz5rlrp5Y1eLC7HqqIiEgYip4WNbiWs8/nLh4+aBCcc06oSyQiInJc0RXUAC+9BPff7+6P1R0uIiISJmql69sYM9gY860xZrUxZnxtfEa1de4MkycrpEVEJCIEPaiNMTHAE8AFQGfgcmNM52B/joiISDSojRb1qcBqa+1aa+0BYApwYS18joiIiOfVRlBnAT+WeZ7rX1aOMWacMSbHGJOTn59fC8UQERGJfLUR1KaCZUfN32mtfdpam22tzc7IyKiFYoiIiES+2gjqXKBFmefNgQ218DkiIiKeVxtB/QXQ3hjT2hgTD4wGZtfC54iIiHhe0M9RstYWGWN+BbwHxADPW2uXBftzREREokGtnExsrX0beLs2ti0iIhJNomeubxERkQikoBYREQljCmoREZEwpqAWEREJYwpqERGRMKagFhERCWMKahERkTCmoBYREQljCmoREZEwpqAWEREJY8bao65AWfeFMCYf+CGIm2wMbAni9iKd9kd52h+HaV+Up/1RnvbHYcHeFydaawO6xnNYBHWwGWNyrLXZoS5HuND+KE/74zDti/K0P8rT/jgslPtCXd8iIiJhTEEtIiISxrwa1E+HugBhRvujPO2Pw7QvytP+KE/747CQ7QtPHqMWERHxCq+2qEVERDxBQS0iIhLGPBfUxpjBxphvjTGrjTHjQ12eumCM+d4Ys8QYs8gYk+Nflm6Med8Ys8p/7/MvN8aYx/z7Z7ExpmdoS19zxpjnjTF5xpilZZZV+fsbY8b6119ljBkbiu8SDMfYHxONMev9vyOLjDFDyrx2h39/fGuMOb/M8oj/WzLGtDDGfGiM+f/tnE1oXFUUx3+HmKRii0n9IrSCiXRhF6UGkYLShUq03UShi6wUFQQ/wC5cFAqiSwXdiQVRqCK2WhW7ES1+4MpU1KamBOu0CpaGZlFbdePn38U9rw7DzEjKm7zx3fODxzvv3Mtwz3/OnTNz700WzOyYmT3u/izzo4seuebHKjM7bGZzrsfT7h83s1l/r/eb2ZD7h/254e3XNb1WW51KQVJtLmAAOAFMAEPAHLCx6nGtQNw/AFe2+J4Fdrm9C3jG7e3A+4ABW4DZqsdfQvxbgUlg/mLjB9YCJ/0+6vZo1bGVqMdTwBNt+m70eTIMjPv8GajLXALGgEm31wDHPeYs86OLHrnmhwGr3R4EZv19fxOYcf8e4GG3HwH2uD0D7O+mU1njrNsv6puBhqSTkn4H9gHTFY+pKqaBvW7vBe5u8r+qxOfAiJmNVTHAspD0GXC2xb3c+O8EDkk6K+kn4BBwV+9HXz4d9OjENLBP0m+SvgcapHlUi7kkaVHSV27/AiwA68g0P7ro0Ym654ck/eqPg34JuA044P7W/Cjy5gBwu5kZnXUqhboV6nXAj03Pp+iehHVBwIdm9qWZPeS+ayQtQpqcwNXuz0Wj5cafgy6P+XLuK8VSLxnp4cuUN5J+NWWfHy16QKb5YWYDZnYEWCJ9ATsBnJP0p3dpju1C3N5+HriCHutRt0JtbXw5/P3ZLZImgW3Ao2a2tUvfXDUq6BR/3XV5Ebge2AwsAs+5Pws9zGw18DawU9LP3bq28eWgR7b5IekvSZuB9aRfwTe06+b3SvSoW6E+BVzb9LweOF3RWFYMSaf9vgS8S0q2M8WStt+XvHsuGi03/lrrIumMfyD9DbzEv8tytdfDzAZJRel1Se+4O9v8aKdHzvlRIOkc8Clpj3rEzC7xpubYLsTt7ZeTtpl6qkfdCvUXwAY/sTdE2uw/WPGYeoqZXWZmawobmALmSXEXJ1PvA95z+yBwr59u3QKcL5YAa8Zy4/8AmDKzUV/2m3JfLWg5h3APKUcg6THjp1nHgQ3AYWoyl3z/8GVgQdLzTU1Z5kcnPTLOj6vMbMTtS4E7SPv2nwA7vFtrfhR5swP4WOk0WSedyqGq03a9ukinNo+T9hl2Vz2eFYh3gnTacA44VsRM2jf5CPjO72vdb8ALrs83wE1Vx1CCBm+Qluv+IH2zffBi4gceIB0CaQD3Vx1XyXq85vEeJX2ojDX13+16fAtsa/L/7+cScCtpCfIocMSv7bnmRxc9cs2PTcDXHvc88KT7J0iFtgG8BQy7f5U/N7x94r90KuOKfyEaBEEQBH1M3Za+gyAIgqBWRKEOgiAIgj4mCnUQBEEQ9DFRqIMgCIKgj4lCHQRBEAR9TBTqIAiCIOhjolAHQRAEQR/zD0hdeV8jV+C3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "from sklearn.preprocessing import normalize\n",
    "plt.figure(figsize=(8,8))\n",
    "np.random.seed(42)\n",
    "\n",
    "\n",
    "class linear_regression:\n",
    "    def __init__(self,fit_intercept=True,valid_freq=100):\n",
    "        '''\n",
    "        Class initilization\n",
    "        \n",
    "        Paramaters:\n",
    "            fit_intercept:boolean,optional,default True\n",
    "                wherether to calculate the intercept for this model. If set to False, no intercept will be used\n",
    "                in calculations (data is expected to be already centered)\n",
    "            valid_freq:int,optional,default 100\n",
    "                How many iterations per validate\n",
    "        '''\n",
    "        self.fit_intercept=fit_intercept\n",
    "        self.valid_freq=valid_freq\n",
    "    \n",
    "    def fit(self,X,y,X_val=None,y_val=None,max_iters=1000,batch_size=32,learning_rate=0.01):\n",
    "        '''\n",
    "        Fit linear model.\n",
    "        \n",
    "        Parameters:\n",
    "            X:array-like, shape(n_samples,n_features)\n",
    "                Training data\n",
    "            y:array-like, shape(n_samples,)\n",
    "                Target values\n",
    "            X_val:array_like, shape(n_samples,n_features),optional\n",
    "                Validate data\n",
    "            y_val:array_like, shape(n_samples,),optional\n",
    "                Validate target values\n",
    "            max_iters: int,optional, default 1000\n",
    "                Maximum iterations\n",
    "            batch_size: int, optional, default 32\n",
    "                batch_size\n",
    "            learning_rate:float,optional,default 0.01\n",
    "                learning rate\n",
    "        \n",
    "        Returns:\n",
    "            self:\n",
    "                return an instance of self.\n",
    "        '''\n",
    "        with tf.Graph().as_default():\n",
    "            dim=X.shape[1]\n",
    "            weight=tf.get_variable(\"weight\",shape=[dim],initializer=tf.random_normal_initializer())\n",
    "            bias=tf.get_variable(\"bias\",shape=[],initializer=tf.constant_initializer(0.0))\n",
    "\n",
    "            X=X.astype(np.float32)\n",
    "            y=y.astype(np.float32)\n",
    "            train_dataset=tf.data.Dataset.from_tensor_slices((X,y))\n",
    "            train_dataset=train_dataset.batch(batch_size).repeat()\n",
    "\n",
    "            val_dataset=train_dataset   # val_dataset must be exist because the tf.cond, see the tensorflow page\n",
    "                                        # for detail. This is unnatural.\n",
    "            valid_flag=X_val is not None and y_val is not None  # Wherether to perform validate\n",
    "            if valid_flag:\n",
    "                X_val=X_val.astype(np.float32)\n",
    "                y_val=y_val.astype(np.float32)\n",
    "                val_dataset=tf.data.Dataset.from_tensor_slices((X_val,y_val))\n",
    "                val_dataset=val_dataset.batch(batch_size).repeat()\n",
    "\n",
    "            training=tf.placeholder(dtype=tf.bool)  # training or validate\n",
    "            x,y_true=tf.cond(training,lambda:train_dataset.make_one_shot_iterator().get_next(),\n",
    "                            lambda:val_dataset.make_one_shot_iterator().get_next())\n",
    "\n",
    "            bias=tf.cond(tf.cast(self.fit_intercept,tf.bool),lambda:bias,lambda:0.0)\n",
    "            y_pred=tf.reduce_sum(tf.multiply(x,weight),axis=1)+bias\n",
    "            loss=tf.losses.mean_squared_error(y_true,y_pred)\n",
    "            train_op=tf.train.GradientDescentOptimizer(learning_rate).minimize(loss)\n",
    "\n",
    "            train_losses=[]\n",
    "            valid_losses=[]\n",
    "            with tf.Session() as sess:\n",
    "                sess.run(tf.global_variables_initializer())\n",
    "                for i in range(max_iters):\n",
    "                    _,loss_=sess.run([train_op,loss],feed_dict={training:True})\n",
    "                    train_losses.append(np.array([i,loss_]))\n",
    "                    if valid_flag and i%self.valid_freq==0:\n",
    "                        loss_=sess.run(loss,feed_dict={training:False})\n",
    "                        valid_losses.append(np.array([i,loss_]))\n",
    "                train_losses=np.array(train_losses)\n",
    "                valid_losses=np.array(valid_losses)\n",
    "                plt.plot(train_losses[:,0],train_losses[:,1],'r')\n",
    "                if valid_flag:\n",
    "                    plt.plot(valid_losses[:,0],valid_losses[:,1],'g')\n",
    "                    plt.legend([\"training loss\",\"validate loss\"],fontsize=15)\n",
    "                else:\n",
    "                    plt.legend([\"training loss\"],fontsize=15)\n",
    "\n",
    "features,targets=load_boston(return_X_y=True)\n",
    "features=features.astype(np.float32)\n",
    "targets=targets.astype(np.float32)\n",
    "features=normalize(features,axis=0)\n",
    "\n",
    "X_train,X_valid,y_train,y_valid=train_test_split(features,targets,test_size=0.3)\n",
    "\n",
    "ols=linear_regression()\n",
    "ols.fit(X_train,y_train,X_valid,y_valid,learning_rate=0.01,max_iters=3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Robust linear regression\n",
    "It is very common to model the noise in regression models using a Gaussian distribution with zero mean and constant variance, $\\varepsilon_i \\sim \\mathcal{N}(0,\\sigma^2)$. But Gaussian distribution is very sensitive to outliers in training data.One way to achieve robustness to outliers is to replace the Gaussian distribution with another distribution which has heavy tails. One possible is to use the Laplace distribution as following\n",
    "\\begin{align*}\n",
    "p(y|\\vec{x},\\theta) &=Lap(y|\\vec{w}^T\\vec{x},b)\\\\ \n",
    "                    &\\propto exp \\left[-\\frac{|y-\\vec{w}^T\\vec{x}|}{b}\\right]\n",
    "\\end{align*}\n",
    "Unfortunately, the obejective function of Laplace distribution is hard to optimize.\n",
    "\n",
    "An alternative to using a Laplace likelihood is to minimize the **Huber loss** function  defined as follows\n",
    "$$\n",
    "L_H(r,\\delta)=\n",
    "\\begin{cases}\n",
    "\\frac{r^2}{2} & |r| \\leq \\delta \\\\\n",
    "\\delta|r|-\\frac{\\delta^2}{2} & |r| \\geq \\delta\n",
    "\\end{cases}\n",
    "$$\n",
    "This is equaivalent to Gaussian loss for errors that are smaller than $\\delta$ and Laplace loss for larger errors.This is  everywhere diffenrentiable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f77d0e218d0>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAHVCAYAAAAgiIjxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xd4VNXWwOHfmZIOBEKTAEEExIA0o6ggTVSUogiiAiI2EBWwAEEBaaJU0SvWq4h+liuCGKpKFURAAZNQlSIhCSUkJKRNpp7vj0kmM8kkmYT0rPd5eC6zc+acPVxhzS5rL0VVVYQQQghRsTQV3QEhhBBCSEAWQgghKgUJyEIIIUQlIAFZCCGEqAQkIAshhBCVgARkIYQQohKQgCyEEEJUAhKQhRBCiEpAArIQQghRCejK82H169dXW7RoUZ6PFEIIISrMgQMHElVVbeDJteUakFu0aMH+/fvL85FCCCFEhVEUJcbTa2XKWgghhKgEJCALIYQQlYAEZCGEEKISkIAshBBCVAISkIUQQohKoFx3WQshhChdqampJCQkYDabK7orNY5er6dhw4bUrl27VO4nAVkIIaqo1NRULl68SHBwML6+viiKUtFdqjFUVcVgMBAfHw9QKkFZpqyFEKKKSkhIIDg4GD8/PwnG5UxRFPz8/AgODiYhIaFU7ikBWQghqiiz2Yyvr29Fd6NG8/X1LbXlAgnIQghRhcnIuGKV5p+/BGQhhBCiEpCALIQQQlQCEpCFEEJUmNGjRxMWFgaAyWRi8uTJ3HHHHTVy17gEZCGEEJVCZmYmn376KX5+ftx+++0V3Z1yJwFZCCFEpRAYGMjly5f5+eefGTx4cEV3p9xJQBZCiDKSlG4kKjaFpHRjRXelyqhp09TO5KQuIYQoAxGR8YSvjkav0WC22Vg4pAODOgVXdLdEJSYjZCGEKGVJ6UbCV0eTZbaRZrSQZbYxZXV0+YyUFaXif4kSkYAshBClLC7ZgF7j+s+rXqMhLtlQQT0SVYEEZCGEKGVN6/pittlc2sw2G03ryjGXomASkIUQopQFBXizcEgHfPQaannr8NZpeL5Xq/J5uKpW/C9RIhKQhRCiDAzqFMzu8D4806MloPLJztN0W7CNtZHxFd01UUnJLmshhCglSelG4pINNK3rS1CANwAf7DiJ0aJitFgAmLI6mm6t6jt+Llxt2rSJjIwMIiMjAVi1ahUAN998MyEhIRXZtTInAVkIIUqBuzSnkCB/9BoNWeSuJ+ds7pKA7N64ceOIiYlxvH7ooYcA+Pzzzxk9enQF9ap8SEAWQoir5JzmlBN8p6yOZv0L3WVzVxFWrFjh8vrMmTMV0o/KQNaQhRDiKhWU5pRhsrJwSAe8dQp+ei3eOoWFQzrI6Fi4JQFZCCGuUmFpTvY9xwoq9g3IaUZLBfRQVAUSkIUQ4irlTXPy0WtYOKQDAOGrozFabBjMVkxWlWlrDvP13pgi7ihqIllDFkKIUjCoUzDdWtV32WUdFZuCTpP/KMlZ6w7TrJ4v7ZrUkelr4SABWQghSklQgLdLgG1a1xeTNf9BGWYrPPvVQWyqKkUnhINMWQshRBkJCvBm5sBQtz/LNFnLt+iEqPQkIAshRBka0TWEeQ+0x0ur4KvL/0+uFJ0QOWTKWghR47k7Yas0jbg1hH7tG3PkXCrPfLkfoyV3R7bkJYscEpCFEDVazglbOo2Cyaoyc2AoI7qW/hGNQQHe9GjTgEVDOzAlz4lesrFLgARkIUQN5nzCVo5paw6Dah/Vlvaz4pINdGtVn93hfcp0RC6qJllDFkLUWHHJBrdpSbPXHSnVjVYRkfF0W7CNkZ/uo9uCbew+mUjHZoESjIHRo0cTFhYGwJ9//skTTzxBq1at8PPz4/rrr2f27NlkZWVVcC/Lh4yQhRA1VkFpSXpt6RWAKOica6n4lN93333HqVOnCA8Pp3Xr1kRHRzNjxgyio6NZvXp1RXevzElAFkLUWDlpSdPWHHZpt6pqkRutPN0IlnPOdVEVny5c+AIfn+sIDOxewk9T9YWHh9OgQQPH6169euHj48PYsWOJiYmR8otCCFGdjegaAqp9mlqv1WDNPqyjsCDrrtRiQYd72EfhVpe2vDurL1z4kuPHn0Cj8aNDhw0EBvYsnQ9XxTgH4xydO3cGICEhQQKyEEJUdzlpSZ6MeIs7Bf3byURsTrPiOg0uAf/ixW84fvwJQMVmy+Dff6fTqdNOFCX/2nZN9Pvvv6PRaLj++usruitlTgKyEEKQ/9jLgng6BQ25wdvstE6t1Wjo1qo+AKfOfk3s6VGQfS9//w60b//jVQVjZXbFB3J1Zv51+ZK4cOEC8+bN47HHHqN27dqlcs/KTAKyEEIUQ2GlFvOuK7sL3kaLjW/2naVlrR34pD+LVmP/mVVzPR07bkGvDyrXz1NZmUwmhg0bRkBAAEuXLq3o7pQLSXsSQohiKKjU4m8nE11Sm9ZGxmevH9vy3WPnoS/xSR/nCMbn0psy9dfXSTNV/1GgJ1RVZdSoURw5coSNGzdSt27diu5SuZARshCiRivJsZl5Sy0CdFuwLd+68u7wPrzQuxVLNv/jeG/HBn8ytsN8tBr7Rq/zGcEs+PNNbASVSqpVaU0XV6SXXnqJiIgINm/eTNu2bSu6O+VGArIQosYqzm7pvJzXnKNiUwpcVx7etTnLtp/AaFFpX/8AL3Seh05jAeBixjUs+GMeV4z10KT54GuVM63feust3nvvPVauXEn37jUrBUymrIUQ1UJSupGo2BSPT9hy3i2dZrRcVSnEwtaVgwK8WTS0I50aRjGx8xvos4OxTdOcd/6aj5XGaFJrkfbNLQy7V+Xy5WI/vtr45ptveO211xg1ahTBwcHs3bvX8evSpUsV3b0yJyNkIUSVV5yR7smLaUTGpuCj1+Q7NrOg3dJFyVlXLqhoxB0h/1An7A1UmxkAb+8QOnf+lfVdGrPzDyPjhnqTmOJNYgrcf4+BnX/4UhOznn755RcAVqxYwYoVK1x+9vnnnzN69Ojy71Q5koAshKjSPM0LTko3Mv3Hw2w6fKHAe11NKcScdeUj564ACu2a2DdopaTs4tChAag2e81jb+9mdOq0DR+fEE4dgXFD4WKyvZ9+ZPBG2ssofATUjIjsHHjdBeKaRAKyEKJK8yQvOCIynsnfR7k9t9rZjAGhV7Wp6reTiS4j9SX3G6mVORKbLRMAL68mdOy4DV/flhw6BHfeYeTSFfvz/ElnU9Bj3BExnxo5PBayhiyEqNoKW7+F3BF0UcHY30tL+yZ1StyPvGvSwX5H0aQMx2bLAMDLqzGdOm3Hz68VUVHQq1tuMK5FKj83eIw79i6CGnAilXBPArIQokorKC84Z6Qbl2xA68GI07mgRHE3iOU8R6+x/5N6be1/eCVsJr46+zS1Xt+Qjh234efXhoMHocftWVxOs/evNldYW28k3fYugVativXZRfUiU9ZCiCovb16w87Tz4fgrZJisbt+nUcDfS+eyCaukqVA5I/WQ2ieZdPMM/PT2kbFGG0SnTtvw97+B/fvhrp4mUjN9AKhDCl/6P8rs4Q9zY8Ng5Iyumk0CshCiWnB3FnVSupG5G47mu3ZI52DG9bqOuv5eLkH8amoXBwV4M/s+8E+bgX92ME4312LpvjlMrBVIIyPc08fElUwvAOpymRW1HuXNkcO54F+fb/adZfydrUvjj0JUUTJlLYSotpynkXP4e2sZdXsLWjWqRVCANx2bBbpMb+e9PmeDWFHWHdiKX9ojBHilAZBuCmDhH29wKjmE8e+c5a5eucG4HkmsqD2cNx4bwYXa9kITy7afKFEOtKg+JCALIaqNvGu/7jZ8WW1qgalNRW0QK0j8pUjUpKHU8koFINPsz6L9b3A27Tqy4upy4avOpGXZg3F9LvG/hqOY+9hIEmrlTlJ7abUeBX5RfUlAFkJUCxGR8fmKOxS14Suv4l4PkJn5NyeP30NtrxT7a7Mfi/bPISa1FVmx9bj8v5swWe1rxg1IYHOrscwYPpJLAa4FE0zWkudAi+pB1pCFEFVeYWu/hW34cqc412dmniQysg+qNQEAg8WXtw/M5t8r15MVE8TllZ0x2+zvb8QFPqk/msyv3iFjYwwYLS73eqF3q6suLCGqNhkhCyGqvKLWfvOuFbvjPN3tyfUGw2mionpjMp0DwIYvSw/M5GTKDRjOBHF5ZRdHML6Gc3zYYDRzH3sCS92gfNPi3jqF4V2bl+izV3WjR48mLCzsqu8za9Ys6tevXwo9qjgSkIUQVV5J135zuJvuLozBcIbIyN4YjXEAaDS+dO64kcd7PYI5piGXv++C2WZfMw4mjg8aPsHrI57kklcATer4MCysqcv9Hr65mYyOhQRkIUTVV9Tab2EHfRS36lNW1lmiovpgNJ4FQKPxoX37tdSt24u6SSEkfd8ZS3Ywbkos7zZ+hunDnybN2x9vnYZzVwys3B/ncs+V++Nkh3UlZrVaMZlMZf6cIgOyoijNFEXZrijKMUVRjiiKMjG7vZ6iKJsVRTmR/b91i7qXEEKUlUGdgtkd3oevnu7K7vA+jsM8ihr9upvu1ihKdpEIV0ZjPJGRfcjK+je7xYvm131PvXp9Wb9OZfAgC0arfWtOc2JY2mQMMx59inRvPwDsh3cqJU6tqs42b95Mhw4d8Pf3p3v37hw5cgSAM2fOoCgK69evd7m+oKnu3bt306VLF3x8fOjUqRO//fZbvms+/fRT2rVrh7e3NyEhISxcuNDtvX/88UfatWuHj48P+/btK8VP654nI2QL8IqqqjcAtwLPK4oSCkwFtqqq2hrYmv1aCCEqTN61X09Gv/5eWoxW1+nuTJOVZ77c7xK8L1w+w779PcnKOgWAxabjw6jp9PtIy7SliTz4gBVTdjBuwb8sDh7H9EeeJtMrd9p85sBQ2jWpfVXT69XR2bNnmTx5MtOmTePbb78lISGBYcOGoaqFnz+eV2ZmJiNHjuTZZ5/l+++/JzAwkHvvvZcLF3IrfC1atIhx48bxwAMPsH79esaNG8eMGTNYtmyZy73OnDnDlClTePXVV9m4cSPXXnttqXzWwhS5y1pV1fPA+ezfpymKcgwIBu4HemVf9gWwAwgvk14KIUQJFFUJKueYTMXNP/xGi8qU1dGEXlObXw4fpm7Ww1zjHwuAxaZl2V+vEnmpC5l/N2BBRCBW1f7PaUtOMb/5BKYPfRqjPndd+LV72zKiawiAS+1kk9XG871K7wzrylAoqphxlMuXL7N7925at7afVGaz2Rg8eDB///03Pj4+Ht/HYDAwb948hg8fDkDv3r1p3rw577zzDvPnzyc1NZXZs2czffp0Zs6cCcBdd91FZmYmb7zxBuPGjUOr1QKQlJTEli1b6NSpU/E+zFUo1hqyoigtgM7APqBRdrDOCdoNC3jPGEVR9iuKsv/SpUtX11shRJVTkkINpaWwzV7Oo2djIZWgHvpgHbUzH3EEY6tNw4dR4URe6krGscYk/djZEYxbcYJf+8zBd+V/UPx88ffW4qXTMG9we8b0vM5xz5zp9Wd6tARUPtl52qPNZNVVixYtHMEYIDQ0FIC4uLiC3lKgwYMHO34fEBDAXXfdxR9//AHAnj17yMjI4KGHHsJisTh+9enTh4sXL7o8Lzg4uFyDMRQjD1lRlABgNfCiqqqpiodfw1RV/QT4BCAsLKyY35uEEFVZSQs1lJaczV5T8vQhKMCbqNiUfKPnvHQk83LYNJrWsm/gsto0fBQ1hQMXbyfjaGMur+uEDfuIqg1/s/2eBTRZ+1+aennRtkVDfjuZSP0Ab267zn3ZiA92nMRoUTFa7DnJnp6bXd0EBga6vPbysm+Ky8rKKtZ9AgIC8PV1nfpv2LAh0dHRACQmJgLQrl07t++PjY0lJMQ+i9GoUaNiPbs0eBSQFUXRYw/GX6uq+kN280VFUa5RVfW8oijXAAll1UkhRNVzNYUaSlNBB324Gz3rtQoaxX6MpVa5wstdZtCs1hkAbKqGTw69wp8Xu5N15BqS1ndEzQ7GbTnGtv5vc82aj0GvJyIynldWRmKx5d53yUMdXb6MFDWdXlLFnS6u7HKmrPPucr58+XK+a9PT0zEYDC5BOSEhgWuuuQaAevXqAbB+/Xq3Afd6p1rUng46S1ORAVmx9+oz4Jiqqm87/Wgt8DgwP/t/I8qkh0KIKqmsAo6nktKNLkE47zPzjp5NVisv9G7Nve0bE598EculIWitpwGwqQr/jX6Jfed7YjocTMKGDqjZK36hHGHbA+/R6PsPQacjKd3IlFVRjmAMYLaqTF7l+mXkanOna4qGDRui1+s5duyYoy09PZ09e/Y4RrPO1qxZ41hDTk9PZ/PmzYwZMwaA2267DV9fX86dO0f//v3L5wMUgycj5G7AY8AhRVEis9tewx6IVyqK8hRwFniobLoohKiKKjLgeDpVnjN6/nrfWd7ffpJPdp7mvzsjebHLdFrW+cdx3fLDE/jjQm+6mTrz7YbGjmB8I9FsGfoxDb9dRlKWlbjkdK4YTGgVDeBag1mrUVy+jBQ2nS5yaTQa7r//fpYuXUpISAiBgYEsWbIk39Q0gK+vL9OmTSM9PZ0mTZqwePFiTCYTEydOBOxT47NmzWLixInExMTQo0cPbDYb//zzD9u3b2fNmjXl/fFceLLL+jegoLH7naXbHSFEdVFRAackU+X2tVwbiprOK2GvuwTjzw+/wG/xd2GMasY3PzVxtHfiLzY/+jn1/+8/RBy64PgCYLJasdryzxu7qzJV3HO2a6ply5YxZswYnnvuOerWrcu0adP4/fffOXz4sMt1fn5+fPnll4wfP55jx47Rtm1bNm7c6JiyBpgyZQpNmjRh6dKlLFmyBB8fH9q0acPDDz9c3h8rH6W4eV5XIywsTN2/f3+5PU8IUfHyTh2XtajYFEZ+uo80p+INtbx1fPV0Vzo2c908dPJiGmujzvHprtNYbRm8fNNMrq931PHzL448x/bY+8jY35TErR0d7V04wOZRX1Hv8yWcvJTBff/Zhclpp7ZOY1/LzWlyt4ZcGo4dO8YNN9xQqvcUxVfY/w+KohxQVdWjw7ql2pMQoky5W78tS55Olb/+4yG+3GvfPe2lzeLlm2a7BOOvjo61B+M/m5G4rYOj/Wb+4Oenvqfuf98mIuock1dFuwRjAF+9jvdHdMl+pdKuSR0Z/YoiSUAWQlQrnkyVn7yY5gjGeo2RF7vMoW293OnPb449zZazA8nc14zEHbnBuCt7+enZCAI/WEhShonw1dGYLPnTpsw2G+2a1JYgLIpFArIQotopam02MjYFAL3GxIQubxAaFO342ZqTT/JLzANk7GlO4s4bHe23s5tNE37C/MZMouKucMVgdpvH7KVVZHOWKBEJyEKIaqmwqfJOzQLRKWZe6PQmN9b/y9HuFTiNn2O6kbG7KYm/5Qbj7uxi4yvb2DZiLOELtzuOvLTmmRr30mnYOL47rRrVKpsPJao1Kb8ohKhxWjbw5q0+79CxYe4m09NZz3F7pze4+eLNLsG4B7+y6dVdmGaGE/7DIUehCqPFhqIoeOsUR8nHxUM7OIKx85GhFXl8qKg6ZIQshKj2nHd61/XTcPToIzTQ/+r4uclnIveHLWD2a0a++m99R3tvtrHu9f34z36NU26O2vTRaXl/RGfq+Hq5TI0750EbzBYURcFHp62Q40NF1SEBWQhRrTkHR6tq5t17PkZvzq2t+9OZh1h/uh+XXzrF+b2hjva+bCbijcP4TZsCFLx7O+8Oand50KBittbs86pF0WTKWghR6eVM+Z68mFasqV/n4JhuNPJY28UuwXhzzIN8e2wUCZuauwTje/iJtQuO4zftJUdbzu5tH73GMUXtbvNWzpGhBck5PlSIvGSELISo1HJGuKpNxWhV8dHbg50nU785wdGImadvfIfbmuROU3vXGUvEqQcwbAnh0sHcQx3uYwOrl8bi8+L4fPfz5GQtdyNpZ3JetSiIjJCFEJWWu5rFWWYbWWYbU1ZHFzlSblrXF4vNwhPt36Nb8HZHe70GY2l13Ttc3tDCJRgPZC1fLjmLz4vPFnjPoABvOjYLLHDKOe9IWqexn9RV2KhaCJARshCiEnNXMSqHJ5Wj6vnreffer/EybXG0ZWgfpU7Am0wbb+L8wbaO9vtZw/OTThL30DOQbryqoJl3JJ3zWeS8aleelDjcvn07vXr1KvvOVAISkIUQlVZh079FTf2qqsqJEy/gZfrW0ZahHcYrP48gbV0i56LbONofZDUhA/7gGW0PvD/6HUWjeLwbuqCzuvPmQUsgzm/Pnj2O3xsMBvr06cP06dNdSiOGhoa6e2u1JAFZCFFpOR+D6W4NuaAgp6oqJ0++yLlzHzraAoNG8ux3j5C8phUJR3OD8TC+I6j/X/zQrgeAfWrcqnq0G/rrvTHMXn8UL62CxaZKSlMx3XrrrY7fp6enA3Dddde5tBckKysLHx+fMutbRZA1ZCFEpTaoUzC7w/uw8tnb2fJSD74bcxu7w/sUGPhUVeXUqVeIj/+Po61hwxHYar1D2po2JBy93tH+CN/QaPBhNra/I999rFaVI+dSC+zX13tjmPbjYUwWG+lGq8fr2qL4PvroIxRF4eDBg9xxxx34+vry3nvv8dNPP6EoCidPnnS5/tZbb2XkyJEubdu3b6d79+74+vpSv359xo0bR2ZmZnl+jCJJQBZCVHo5G6laNapV6IYqVVU5fTqcuLiljrYGDR6mTZsVLJ6o5dzR1o72EXxFnSHH+Tn0Nrf3MttUnv7iT9ZGxuf7WVK6kdnrjuRr12oUSWkqQw8//DBDhgxh48aN3H333R6/b9u2bdx99920aNGCH374gcWLF7NmzRrGjBlThr0tPpmyFkJUC6qq8u+/04mNXeRoq19/CNdf/xVPDMnk27W5tZAfV1bgNTSG+14byx1ZFqb9eNjdLTEVMHUdl2xAr9Vgslpdrjdb1QpPadqxo+iNUmWtVy+16ItKYNKkSYwdO9bx+vz58x69Lzw8nL59+/LVV1852ho2bMjAgQOZOXMmrVu3LuTd5UcCshCiWjhzZhZnz77peB0UdD9t2nzLqAcy+XZDbUf7U9rlTF1eh7oPTiUowJuo2BT8vTRkmNxvHnO3m7tpXV+sav6gM3NgqGzeKkPOm708lZKSwoEDB/jss8+wWCyO9p49ewJw8ODBShOQZcpaCFHlnTkzl5iYOY7XQUEDaNNmJSMHuQbjsdpP+WRTCK1GDXEETntwLfje7nZzO+ca+3tp8dIqzHugPSO6hpTuBxMuGjVqVOz3JCUloaoqTz75JHq93vErICAAm81GbGxsGfS0ZGSELISo0mJi5nPmzOuO1/Xq9aN16+8ZMSCT1Vtyp6mf133Me5vbovTq6fJ+553cWkUhy2IFFfy8dI5iEO5GvZ6c2lURymq6uDLIm7ecs8vaZDK5tF++fNnx+7p16wLw1ltv0bdv33z3bNq0aWl3s8QkIAshylRBebqlITZ2Cf/++6rjdd26d9G69Q88cm8WETtyg/FE/Qcs3dYRpXs3t/fx9CCPvJ+lsJrLouzlBNNjx4458pVPnTrF6dOnueWWWwCoV68enTt35sSJE0ydOrXC+uoJCchCiDLjXGmptEsPxsW9y6lTkxyvAwN707r1jwy728j633KD8cte77F4x80otxWe21rUQR5l+VlEybRq1Yobb7yRV199FZ1Oh8lk4s033yQoKMjlukWLFnHvvfdis9l48MEH8ff358yZM6xfv56lS5cSElI5lhpkDVkIUSacz6FOM1o8ytPNqepUVC5vfPz7nDz5ouO1f0B3Wrdex9C+JpdgPMX7XRbvurXIYFzU80vyWUT5+O6772jUqBHDhw9n5syZzJs3j2uvvdblmjvvvJPt27cTFxfHiBEjGDRoEEuWLKFly5b5gndFkhGyEKJMuDuHurDzpz0dgZ479wknTrzgeP1PcihLf3qJWuvT+fNI7qafV32WMm93L5QunT3qb97nzxgQSvsmdWha17fYn0UUX0BAAKqbnesAzz77LM8+677gxw033MCuXbtc2u67775813Xr1o3NmzdffUfLkARkIUSZcHcOdUHnTzuPQHOCnrv83/Pnl/PPP7l5qCdTrmfx3jmkftGZ4xdzg/EM38XM3nM3SscOHvXV3fOnrTmMv5cWq6oyY0Cox59FiJKSKWshRJnIW4awsNKDOSNQZzkj0BwXLnzB338/7Xh9+kprFv0+jysrbuLixdwpytn+C5nz570eB+OCng+QYbIfiTl3/VFm9A/16LMIUVIyQhZClBlPU4OKGk1fvPgNx48/AdinNGNSW7Hw93mkrAgj4VLuhpy5AfOZ/ucD0LYtxVFYVSmwfzloH1yH3eF9Kl2ak6g+ZIQshChTOedQFxbAChpNJ2eYWLvnfY4de4ycYOzv3xGfemtIWd7VJRjPCXiT6QeHFDsY532+v7c2389zvhx48lmEKCkZIQshypS7PGR3bXlH0+9u+Ycvt7/PuI4LQGMfvWq8bqBuow28e7cvFxPrOZ4x3XsWt317N1zFEYjOzz8cf4W5G466bDCrrEFYVdV8B2aI8lPQRrSSkIAshCiQc+AE94dlFMbdzmkV3O6mdn5WcoaJQ6dX8nynhWizg3F8WnPe2TeL2A9sXEzODcav+8xgy+iWfL/vCrt7Ga8qcObkIndsFki/9o0r/fS0Xq/HYDDg5+dX0V2psQwGA3q9vlTuJQFZCOGWczDNslhRVRVfvc7jQzHc7VyevCoaUDFaVJfd1GlZFuZuOIoGMFltPHHTCZ7vNB+dxl5N6Xx6U97atYCET27n0pXcow5n+U7jp9GtOV+7AbVKOQ2pKpzC1bBhQ+Lj4wkODsbX11dGyuVIVVUMBgPx8fElOmPbHQnIQoh83AVTgDSjvVqOu5SkvNzl7mo1CqgKYHVpm73uCKbsCg831j/ALXXnotPYn3Uhowlv7lpIwsfduZSaG4zn+r/K2sdvIKGW/WCHmpiGVLu2vXDGuXPnMJvNFdybmkev19OoUSPH/w9XSwKyECIfd8HUmSeHYrjbuWy1qeRszsphtqposwd27YL+YkLnN9BnB+OEzMa8uXMRFz/S0jG8AAAgAElEQVTswaU0+4hcwcbcgFdp+91QvtqTSK0qsM5blmrXrl1qAUFULAnIQoh8ikoD8mQ06lxFyXm9GHBpmzEglJkRh7mhXhQTu8xFr7WP9BINDVmwczEZy/tyKa0hYA/Gb9QKp8uPI+nXpyO9epRd4QohypsEZCFEPnmDqbs1ZE8CYEF5yHnbbIbdNLTMwUtrL6OXZGjAmzsWEfdBHy6l5wbjT4Jn8OBvL1KvRbCjn3n7kZRu5Mi5VEClXZM6EqhFlaGU5pbtooSFhan79+8vt+cJIYqvoJ3Vzr8vzSCXkvIb0dH9sNkyALicFcQb298mdtldJGZeA4AGK1+2mMmIAy9DvXoF3isiMp5XVkZiyR7c67UKSx7qKFWZRIVRFOWAqqphnlwrI2QhhENRBR48DcSe1kC+cmUPhw7d6wjGFhqycOdCYpfdTWJmYwC0WPiq1Swe+XMSBAYWeK+kdCNTVkU5gjHY16cnryp6A5oQlYEEZCEE4HmBh7zvyRt4IyLjmbIqGq1GwWpTWTTUfYpUauofREf3w2pNB0Cvb0Tza3aQ+HEDEjPtO6d1mPm/1rN4ZH84FLFxKS7ZgFbR4LyDG+y7uKUqk6gKJCALIYDSKZfYrVV9Jn0fhdmauxT2yvdR+YJ6WtpBoqPvwWpNtT9H34BGjbZyV1hDTmcf+qHDzLwGk1k4uAd+p9MY1KnwgNy0ri9WNf9GNKtNrXHpUKJqkrOshaimktKNRMWmkJRu9Oj6kpZLTDNayDLbmLwqij2nklyCMdinjY+cu+J4nZYWSVRUXyyWFAB0uiAa1N/CXWFNHMFYj4m5DSbx+cg7uKx4M2V1dJGfIyjAm0VDO6Jz+ldNr1VYNLRmpkOJqkdGyEJUQ0WtBbtTUJpSYeUSnUfTRovK9uMJBdzdnmicnn4oOxgnA6DT1aV+0GbuuqUZZ1PrAuCFkdkNp/D5yJ4Y9fZnezrtnLOrW3ZZi6pIArIQ1UxJ1oJzFKdcoslqzde+LvocWgWcB8k6DbRrUpuMjKP8FXknVksSAFptHerV/YW7bmlBXLo9GHuTxZxrwvlseC+MOi/HPcxWz6edgwK86dGmgUfXClGZyJS1ENVMzujVWc5asCc8LZf4Qu/8lZW8dVpe7NsGb52Cn16Lt07h7WGd8FH+Zd/+XlgtlwAwWPw4kfwtd3Vt5QjGPhiYFRzOipG9XYIxwMyBoTLSFdWejJCFqGaKsxZ8NYZ3bc6y7ScxOuUZmW02hndtzvCuzR2jbF/NWQ7+1QeNmhOMfZm/ZQHRS7qSYrGnMfmSyfSmU/n0kTtR9Hq8FdBpNJitNmYObMeIriFu+yBEdSIBWYhqpjhrwVf7nEVDC35OUIA3BsMp/vqrNxbzeQCyLD4s2jKfyCUjSLXYR8Z+ZDCt2assf/hOLFodtXRa3h/RmTq+Xo5pc0/zmoWoyiQgC1ENeboWXJbPMRjOEBnZB5MpHgCT1ZvFW97iwOLHSLXag7E/6bx27TQ+G9oXq0YL2EfZzpuxSrJBTYiqSNaQhaimPFkLdsddulRhKVTunpOVdZaoqN4YjWcB0Gh8iE3+kj8WP+EIxrVI5ed736X996+g99ZTy1uHj17jMsp2l17lSQqUEFWRjJCFEA7uRqMquLTN6B9K++A6BY68s7LiiIzsTVbWGQAUxRsf/Q9MGdGdNGstAGpzhZ8HfcCta14FjYZu1zdyO8q2n76luNzfk9KPVVFMSgz74vcxrN2wiu6KqCASkIWoAU5eTCMyNoVOzQJp1aiW22vcpUtNXhUNqBgtqqNt2o+HCfDWYrGp+aaPjcZzREX1ISvrNACK4oW3bhX39LyDFHMAAIEks+6Bj7h1dThk7wZ3V7UJ4HD8FTJMrulVJquNKwYzSenGahGUTVYTS/csZfavs7GpNjo37kzroPw72EX1J1PWQlRzr/94iL5LdzJpVTR9l+7k9YhDbq9zly4FoEHJ15ZutOabPjYaLxAZ2QeD4QQAiqLHS/kf/Xr1cgTjulzmpfZv8HS7jqyNPl9ov5PSjczdcDRfu8Vq4/mvD9JtwTbWRsYX/uEruZ0xO+n8cWembp2KwWLAaDXywqYXKrpbooJIQBaiGjt5MY0v9551aftyz1lOXkzLd627dCmjxYbBkv986Bw508cmUwJRUX0wGP62/0DRcSn5U/rdeTdXLPZgHEQi49rNY/l9vcmyqEWuBRf0BcGqUuXXkxMzE3ky4kl6rujJ0Uu5Xzo6NOrArJ6zKq5jokJJQBaiGouMTfG4PSjAmxn9Q9EX8K+Cn5sfmG02GgdkcOCvPmRmHgNARcuHm2fw+INDSLX4A1CfS4ztuICv+/eG7DVh1aYWeFhJUrqRKwYzJmvBXwageAeeVAY21cZnBz/j+mXX83nk5472AK8A3r77bQ6MOcBtzW6rwB6KiiRryEJUY52aua8f7K49IjKeuRuOotNq8o2U/b21zB7YDqPVxtz1R9FrNJisNib0qk9U9F14qfZgbFM1fLLlNda9NQmD6gdAAxL4dvhXPNW0pyMYAxitKv5eWrf9yNlEZrXZ0GsVfHRaTFYbVpvNtd5xGRx4UlYOXTzEuA3j2B2726X9wRse5N1+79K0dtMK6pmoLGSELEQVVlRFp1aNajHqtuYubaNua55vY5fzhi6DOf+o1GJVaVjbh37tGrM7vA/P9GiJny6VgIzhTsFY4b+bpxLx1lRHMG7EBZ6+5R2SX3oInzzB10evIcNkdfkMedOcLDZQUHl/RGd+n9rHfgynXuM2RaqySjelM/mXyXT+uLNLMG4R2IINwzewethqCcYCkBGyEFWWpwdmzLn/Rkbd2qLQXdbuqjcB+HlpMWePTJ//+qA97WlAKMt3RTGh0wya1z4J2IPx8i1TWDN/OkbVPmJtzHkeufk//HB3D54sYKR+OP4KD3+yx/EZnu/VKl8/TFaIir1CjzYNy+3Ak9KgqioRf0cwYdMEYlNjHe16jZ5Jt09ieo/p+On9KrCHorJRVFUt+qpSEhYWpu7fv7/cnidEdZWUbqTbgm1kOY1mffQadof3KVGQcnc/b53Ckoc68cr3kRgtuf9O1PbO4sUu02lZ57ijbfmWyaycNwcjPgA0IZ7h3T5iXe/bHV8U1kbGuxyzOWNAKHPXH833TFUFU56ayt46Db9PLdlnqwgxKTGM3zSedf+sc2nvGdKTD/t/yA0NbqignonypijKAVVVwzy5VkbIQlRB7ka0JT0wIydH+eW+bXh7yz8uI+5m9fzw0moxWiwAeGsNTOw80yUYf7HlJb6bNxcT9uc2JZYNU7egPvcKU51GsXlHt+4+g5dWy6COTfj6D9ed4V7aqnEYiHNOscGSu9msvl99Ft+1mFEdR6Eo+dPIhAAJyEJUSaVV0en1Hw+5pEUNC2vKiK4hLkUdcp7jpc3ipZtmc11gbprOyu0T+XreW5izg3FzYtg+dzctpz/hKAgBuBSccA6q7j7DE91asOpgXL4qUpV989aumF08u+FZlzQmgDFdxvBW37eo51uvgnomqgrZ1CVEFZRT0elqNji5y1FeuT8Os8VKXLLBcRLWwiEdqOVt5pWb3qBtvcOOa1MS5vLpnIWOYNyCf/l1/l5aTh9ORGQ8t8/fyqOf7OX2+VvdHuBR0Gdo1agWi4Ze3WcrTzk5xT1W9MiXU/z7k7/z8cCPJRgLj8gashBV2NWUJVy1P5ZJq6Lztes04KvXOaat+98YxF9RA0hP3eq4JjNxNvc/9CoW9AC05BTb346k+UtDSEo30vXNLS7pSToN7Hutr9s+FvQZKnvJRZtq4/O/PmfKlilcNlx2tPvr/ZnTew4Tuk5Ap5FJyJpO1pCFqCGcp4CLG8AKylG22OwnYQG89sMBmtjedwnG6QnTeeDh17Bm//PRihNs/89hmo4fAsCRc6nkPdzLYrO392jToNDP4El7ZVBYTvE797xDszrNKqhnoiqTgCxENVCSmsE5Ocpf7smdttYq9qMp7b83M7bDfNJT9zl+nnYhnMGPznQE4zb8zfaP/qHJ2MFOdy5o1q38ZuPKSropnTm/zmHp3qVYbBZHe4vAFiy7dxn92/SvwN6Jqk4CshBVnLsqTVNWR9OtVf0CR5g5o+mJd7Zx5Ci3CPJj5PI/sJptaBULz3VawI31c4Nx6vlXGDx8HjbsB3zcwDG2fvov1zw10OXe7ZrUQa9VMDulLum1Cu2a1Cntj16uIo5HMH7TeMkpFmVGArIQVVxxU6DcjaaHhtmnWBcO6cDUH/7imfaL6dxwr+M9V85NYPCIBajZwbidcoStK+JoNOq+fPcPCvBmyUMdmbwqGq1GwWpTWTS08m7KKsqZlDNM2DRBcopFmZOALEQVV5wUqKJG0wM6NKKZ8iVXkn9zvCc57jmGPPa2IxjfqBxi6zcJNHjkngL7VJVO1CqI5BSL8iYBWYgqLijAmxkDQpm97ih6rX1EWlCaUGGj6Xr+Oo4ff5wryd87fpYc+zQPjloG2TWROyrRbPk+mfpD7ixyE1ll3pRVlF0xuxi3YRxHLh1xaZecYlGWJCALUQUUFvwiIuOzKzApmC02Zg5sV+CGroJG08GBXhw//iQJCd842i+fHc2Qxz8hJxh30fzF5jUZ1BvUs0SbyKqCxMxEpmye4lIaEew5xR/1/0hKI4oyJQeDCFHJRUTG023BNkZ+uo9uC7a5HLLhPAWdYbJisqrM3XC0wOpPbg/jeLA9ifETuHjxS8d1cScfZcjjy3EEY+UA1z70G781vzZfRaYss40pq6MLfGZVUFCdYn+9P0vuXiJ1ikW5kBGyEJVYUWu+cckGtHnWMQva0JUzyu7Wqj67w/sQl2wgONCHpHMvcv78Z47r9vzZn2nhX5ETjG9S/qTxQ3vZ36Ilh1dH88ljYaV2jnZlIDnForIoMiArirIcGAAkqKraPrttFvAMcCn7stdUVd1YVp0UoqYqagf14fgrZJisLu8xmC35NnS5m2Ie2LEJJ06M5/z5jx3X7dvfj2nha1FV++TZLcpegobt53CLlgBoNQqglso52hUt3ZTO7B2zWbp3KVY1989QcopFRfFkynoF0M9N+1JVVTtl/5JgLEQZKGwHdVK6kTnrj+R7T96dv+6nmKM4fGwi586977gu7sQgXgtf7wjGtyq/E/joQY62uDb32VaVdk3qXPU52hUt4ngEoe+HsnjPYkcw1mv0vNb9NY48d0SCsagQRY6QVVXdqShKi7LvihAir5w13yl5RrdBAd5ExaagVTSA6whZn12qEOwj7CsGU55RtsrQ1p+TlLDa8Z5LZwby+LOrsdnsqU23aX7Df/hh/g4Ocbn3zIGhBAV4V9m0JqlTLCqzq1lDfkFRlFHAfuAVVVWT3V2kKMoYYAxA8+bNr+JxQtRMBQW/pnV9saq2fNdbbSqH46/w8Cd70Gs0mKxWbI5Ds1SGtP6Svs2dg/G9PPLUD9hs9n8Oumt2EvjECY43ao5OVfHWaTFb7bu3R3TNDdBVKa3JbDXz9p63mbNzDpnmTEe75BSLysSjak/ZI+T1TmvIjYBE7IfTzgWuUVX1yaLuI9WehChdayPjeXllpKOYg16rMGtQO+auP0qWOTdY6zSg1Wh44Lqvuffarx3tCf/ezfBn1mO12qs29fX6lS+21uVSSHPHmnBVGwXntTNmJ+M2jMtXp/iZLs8wv+98ySkWZarMqz2pqnrR6WH/BdaX5D5CiMIVln+clG4kJMifnyb24NyVLFINZmr76kk1mNDkGe356nUsG7QdNc05GPdhxJh1jmB8j/cO1vwRjG+H1jRxem9VDcSXMi4xZcsUVkSucGmXnGJRWZUoICuKco2qquezXw4GDhd2vRCi+Ao7fCPvz4bd1JSVB+IAXEbGOfo0+xY17QvH64QzPRk5diMWixcA9/lsY/WBFviEtiyHT1a2bKqN5X8tJ3xLuNQpFlWKJ2lP3wK9gPqKosQBM4FeiqJ0wj5lfQYYW4Z9FKLGKSz/GMj3sy/3ni3wXgNarmFwK+dg3I2RY37CbLaPfAf4bmbVX63xvr5FsftY2aazD108xLMbnuX32N9d2iWnWFQFnuyyftRN82du2oQQpaSw/OOc3zv/rCD9W65laJvcv64JMbfy2NjNmM0+ANyj30TqE4n8bAhlUDH6V9mOzpScYlEdyLyNEJVQURWcDGZLkfe4s/k6HmrzieN1wplbGPXsVkwm+z3u0W8g6ckkLtUKcls/uaARcEnqL5elguoUT759MtN6TJM6xaLKkIAsRBkrydRuYfnHSenG7BSdgjMk+jbbyMjQ3BO4EmJu4vFxWzEa7cHpPu91XHgilaQ6QUD+oy8LGwEXt/5yWYlJiWHCTxNY+/dal/YeIT34sP+HhDYILbe+CFEaJCALUYZyAptOo2CyqswcGOqSywsFB+yC8o/jkg346LSYre5HyT2a/szIdh84Xl8624nR47aTlRUAwCN1NnJmdAZJPoGOa5xH30WNgItTf7ksmK1mlu611ymWnGJRnUhAFqKMOAe2HNPWHAYVRtxqD8pFrcW6O3zDXUDM0T14C6PbLXO8vhTbgdHjfsVgqAXAY3XX8/mxW9lw3uh29A1Fj4ALG72XtYLqFEtOsagOJCALUUbikg3oNPlHarPXHaFf+8ZA/t3SnqzF5gTESauiMVlyg+Zt12znyfbvolHsU9mXYtvxxLhfycysDcDooLV8eqw72gb1GNSIAo++9GQEXN5HZ0qdYlETSEAWoow0reuLyZp/ndf5rOnirMU6T20P6hRMkzo+DP14LwBdG//KMx2WOoJxYlxbnnp+JxkZ9mnppxpE8Mnxnmjq5U5TF3T0pacj4PI4OtOm2vj8r8+ZsmWK5BSLak/+SxaijAQFeDNzYKh9mtqJVVUdo01P12LdTW2HBPnjo9dwY72djOmwBI1iv9eluNY8/fwu0tLs07fPNFzDR3/fiSawtsd9rwzFIySnWNQ0npRfFEKU0IiuIcx7oD1eWgV/L61LqcKckWhRZQzdl0+Mxt9LS6cGvzO24yK0GnswPh97LWPG/0Zqqv0AkYcCvmX/KBvrz6QVu+9BAd50bBZY7sE43ZTOlM1T6PxxZ5dg3CKwBesfXc/qYaslGItqSUbIQpSRnCnmfu0b069942LtpHbmbpOVBoX4ixGM7bAAJbv84oW4Fjw/cQ8pKQ0BGBrwDdFPe5Ol9SlWnnBFnsAlOcWiJpOALEQZKM5JVkWtxbrbZNWqzl5MifPQa+ypTwnnmvP8xN9JTm4EwLDaXxH5lC9GL/uJXFabypFzV+jRpmGp9bs0FVSnWHKKRU0iU9ZClLKCpphPXkwjKjaFpHSjx/eJik0BYOGQDnjr7Du22wUdZHznNx3B+NKFYJ6fuJvLl68BYHyz7/jraX9HMAYwW1We/uJP1kbGF7vfnva3JMxWMwt3LyT0g1CXYFzfrz4r7l/Bjsd3SDAWNYaMkIUoZe6mmFWbyn3v/Ya31rORp7uR6n9HhfH2huU81+EN9FozAIkXm/DCxN9JTGwKwKN1v6Dnt9345acTkKfqk8laeFpVeZ/AJTnFQriSEbIQpczdFLPRqmKyeDbyPHkxjcmrXEeqk1ZF423dy3MdZuGlNQFwOaExL0zcTUJCcwAerfc50WPr0rFVowL75lygwpN+l8UJXImZiTwZ8SQ9VvRwCcYdGnVg95O7+WTgJxKMRY0kAVmIUpZ397SXVsFH7/pXraDAGBEZz33/2eVy4AdASMAhkuMfxktrD+LJlxox/qXfuHixBQAjGywnemx93nz4Jlo1qpU9xZ3/r3dhAdbTXd8lZVNtfHbwM65fdr3LAR/+en+W3L2EA2MOcHuz20vlWUJURTJlLUQZcN497e+lZcCy31x+7i4w5qzh5j1M5LrAY7wcNgsfXRYAqZcbMP6lXZw7dx0As9r/jwfWDqJpg1qO4Jnz/G/2nWXZ9hN4abUeHXFZVvnHklMsRNEkIAtRRpx3TxdWuSkn+Llbw21Z529euWkmvjr7aPrK5SDGv7iL+PjWAMzvspLwPx4Crdbt88ff2ZrhXZsXK8CW5glcUqdYCM9JQBaiHLgbeebduDVjQKjLGm6L2ieYFPY6fnp7RaO0lHpMfHknsbHXA7D4lpW8smcoaApfeSqPIy7dKSineNLtk5jeY7rkFAuRhwRkIcqJc2B0V+Jw7vqjzOgfytwNR7m29mle6DQDP30GAGlX6jLx5R3ExNhTgBbeupJXdhcdjCuC5BQLUTISkEWNUx4nURX1jIJSjNoH12HL+AacOj4S1ZYOQEZqHV56ZTv//nsjALNu+orJv4+ASlbzV+oUC3F1JCCLGqU8TqLy5BkFpRjV9/mXM//ci2qzVzbKTK/NS5O2cepURwDe7fU9E7ZVvmBcWE7xW3e+RZBfUAX1TIiqo/LNdwlRRsrjJCpPn+EuxWjxA76c+acfZnMiAIb0Wrw8aQsnTnQB4JN+PzBh29BKFYw9ySmWYCyEZ2SELGqM8jiJqjjPcN7o1cA3jjP/3I3JnABAVqY/r0zZzN9/34yCjU8HruXJiMGFBuPyLApRWJ3i2b1mM/HWiVKnWIhikr8xosYorZOoCgt8xX1GUIA3vppYIiPvwWQ6D4DR4Mfk8J85dqwrCjY+f3A9j69+oNA+OU+Tm6w2XujdiuFdm5dJYD508RDjNoxjd+xul3bJKRbi6siUtagxSuMkqojIeLot2MbIT/fRbcG2fMUaggK8mTEgFC+dBn9vbaHPSEo3cvB0FH9F9sZkOgeAMcuXKVM3cfhwNzRY+XLYBh5fPajQPuWdJjdabCzZ/A+3z99aaDGJ4ko3pTP5l8l0/rizSzBuEdiCdY+ukzrFQlwlGSGLGuVqTqJyl6qUU6wB7NPVh+OvMHfDUfQaBbPFxsyB7dxuGouIjGfBhs280iWcIF/7NLUxy4epr24gOroHWix8+vBGRv6v8GCc89y80+QARotarDrIhZGcYiHKngRkUeOU9KCMgtaHv953lg92nESnUUg3Wl3eM3fDUfq1b+zyvKR0I/PXb+WlLlMdwdhk9Oa1aeuIjOyNFguD23/GW62ak7b1RJFTz+6myZ37dzVr5AXlFPcM6cmH/T/khgY3lOi+Qoj8ZMpaCA+5C3wmq5X3t58ky2zLF4zBfRGJs5dO8VKXV2nodxEAs8mL6TN+5ODBvugw80DH5fzZv6nHU885U/HFLSZRGLPVzILfFhRYp3j749slGAtRyiQgC+Ehd2vQL/RujZe24L9GeQOi0XgBw4UHaOhnXzO2mPW8PvMH/vyzH3pMDOy8gv39XKe4c6aeC0vPGtQpmN+n9uGVu9rgrVPcrpEnpRuJik0pMs1rV8wuOn/cmalbp7oc8PFMl2f4+4W/ebzT43LAhxBlQKashSiGvGvQAO/vOJnvOn8vLVZVdQmIJtNFoqL6YDKeAMBi0TFr9vfs3dsfL4wMCPs/DtzZxO1z9RoNR85doY6vV4Fr34UVk/DksJLEzESmbJ7iUhoR4MaGN/LRgI+kNKIQZUxRVbXoq0pJWFiYun///nJ7nhDlYW1kvEslpxkDQmnfpI5LQDSZEomK6k1GxmEArFYtc+Z8x86dQ/Ami/63fsXRvsEYzTY0GgWLzfXvpV6roFFwKaPo6QljSelGbp+/FaMl954+eg27w/sQFOBdZE7xhK4T0Gv1V/vHJESNpCjKAVVVwzy5VkbIQpSAcy5yUTu3zebLREX1dQRjm1XDG298w86dQ/DBQN/bvuVAj2vAbF+fVlSVcT1bsnz3v3hptZisNqw2G0YrGC0WgGLtnv5631mXYAy5a9vnMv5xm1M8uO1g3u33rqQxCVGOJCALUQxJ6Ua+3neW97efxEvrOv3rLjiazclERd1FRkYUYA/Gb771f+zYMQxfMunT7VsOdW/s8h6rCrddV5+n72hJXLKBKwYzz399kDSjxXGNp7unk9KNvL89/5R6ljWDT6Pn8tHB/2Cx5d63RWAL3rv3PQa0GVCsPxchxNWTgCyEhyIi45myKsox2jRaXHOR84+MU4iOvpv09IMA2GwKCxZ+ztatw/Ejg549vuPwba7BOJfqSM9KSjeW+ISxuGQDXlqNo68AmZq9WAKWs2z/OUeb5BQLUfEkIAvhgZxDQfJO/YL70arFkkp0dD/S0nL3TCxe/Cm//DIKf9Lp3WcVZ7oH4221YrGB1WnNWK9VaNekjuN1zu7uKXk2ZXkyXe2cqmVREris/xiDdh84bbSWnGIhKgcJyEIUISndyPbjCeg07lN98o5WLZY0oqPvJS1tn6NtyZKP2bTpSQJIY/GkvQyd+ShxyQb8vbRsOnyBZdtPotUoWG0qi4bmD7YlPWEsKMCbNwffwNgf55Ck+QZVyY3EUqdYiMpFArIQhchJF9IqChmm/Ad/eOs0zBgQ6jj8I9DXwqFD/UlN/d1xzdKlH7B+/Rhqc4Vb+kaw1KcBQ4EzSRmOVCRQGdez8IIQJTlhbFfMLl7fO45Ebf46xfP7zqeeb71i3U8IUXYkIAtRAOezq5356DXYbDbG92lDPX8v5q4/mh1UM1l61yJ0lj2Oa997713Wrh1HHVIIu3stJzoHUUuj4ci51HznYr+/4yTDuzYvlb5LTrEQVY8EZCEKcORcKhrcTOWqoCga6gXYg3GW2YZVY2Bil7noLJGOyz74YAk//DCBQC7Tqd8GTnYMAshe01XLpDZzTk5x+JZwkgxJjnZ/vT9zes9hQtcJUqdYiEpK/mYK4cbXe2OYte4IZmv+TVxZ2TuWZ687ik4BvcbE+M5v0r5+bjD++OP5fP/9y9QliQ79N/Jve/vUsJfOfpxluyZ1SqU2szOpUyxE1SZnWQuRx9d7Y5j242G3wdiFqmKyGHm+01t0aHDA0fzpp2/wv/+FU48k2g38iTPtc9dpJ93dxpGzfLW1mXNkmDKYsnkKXT7p4hKMQ+qESJ1iIaoQGSEL4SQp3czTBNoAACAASURBVMjs9UfztfvoNFhV1SVIW20mnu80n04N/3S0rVgxk6+/nkZ9LnH9/ZuJbRvocp8lv/zNkC5NCQrwvqrazDkijkcw4acJnL1y1tGm0+iYfPtkySkWooqRgCyEE/tBGgomi2u7VVV5uvu1LN99Bi+tBovVyDMdFtGlYW5q01dfvcYXX8ykAQm0HryVuDZ1yEuruK4Tl7Q2c0xKDBN+msDav9e6tPcI6cGH/T8ktEFose8phKhYEpBFteV83rSnQa9pXd98hR0AbDaVr/aeBVTG9GjOLYHTMaTlpjb973+T+eyzN2isSaDtsB38G1Lb7f2t6tWtE5utZpbuXcrsX2e7lEas71efRXct4vGOUhpRiKpKArKoljwpN+iO86lYWo2CyWLDalOxqpBmtKBgxXDpOQzaHY73fP/9S3z88QKu0VxgxuIzLLgY4PbeOg0sGtqxxLuod8XsYtyGcRy5lD+n+K073yLIL6hE9xVCVA4SkEW1k5RuZMqqaIyW3Bzf4lRHylnb/XrfWd7bdoKcAbOCladufJdbGu9wXPvDD+P54IMlNNWe58c1GYz4IyXf/UbfFkKfGxrRpI4PGSYrSenGYgXlxMxEwjeHszxyuUu75BQLUb1IQBbVjr3coGtKUXFzfJMzTCzbdsKxiUvBxhPtl9E9eJvjmoiIcbz33rs018bz4zojmvYN0O0/5XIffy8tg7s05UxSBmP+b3+hNZPzkjrFQtQsEpBFtWIvN3giX7vJ6vnabURkPJNXRWNyCsajQj+gR9PNjmvWr3+Gd99dRlNNLM1G7efR330ZlNGEdKPr8ZpWVcXfS5vvVK5paw7j76XFqqpup9MLyimWOsVCVF8SkEW1Yt8lrcVocd0m/ULvVh6NjnOOyzQ5RtgqI2/4mN7Nf3Jcs2nTE7z99kc0U85Sd8R+4hr6gtnGyv1x+e43Y0AoGSZrvlO5AMfZ2M7T6RmmDGb/Opule5e61CkOqRPCsvuWSZ1iIaoxCciiWnEuN5jDW6d4fEZ0XLLBKXiqDG/7CXeGbHD8/JdfRrJ48X9pqT1Lg9F/cT6o4FG3v5fWMS2dt0/OcqbTf4v7yW1O8aTbJjGj5wzJKRaimpOTukS14u4ErJydzUnpRqJiU0hKzy1BmLctN3iqPNL2U+5usc5x7datj7JgwQpClBjeXmUkubFPoX2xqqpjjTinT/7e2nzXZdrOM/XXx3jguwdcgvEdze8gcmwkb/V9S4KxEDWAjJBFuStJfnBxuDsBKyIynimrotAqGqyqjUVDO5KWZWH2uiPotRqXtdxhYcEYkubQr0WE4547djzEm29+SQvlNIFPHqHbnfeyMCSAKU6pVcPCmrJyf5xLqlXO53Pu0+H4K8zdcBSdxsZFdTVpXv/j9CmD41mSUyxEzaSoahHn9ZaisLAwdf/+/eX2PFH5lDQ/+GokpRvp+uYWnDdeK0De//J99BrWPd+N/24aw33Xfudo37nzQebM+R8tlDPUeeoIM5+92RFc/b20ZJisjsCf98tGQV8+1h/fyos/v8CplOMufZCcYiGqF0VRDqiqGubJtTJCFuXGub5wSfKDS+rIuVTyZEHlC8YAWkXhn9OzXILx7t2DmDv3W9rqz7DsJx+atr2TjYcvMHnVNry0uV8qOjazn1ntfBSmuy8ft7fxlpxiIYRbEpBFuXHdMGVXGjWAi+bZLFDfZt/ib/rS8XrPnv7Mnr2Slspp1vwRxFGribHv7cJosd8vJ9fZ3ZeKvF8+VGyM+WERJv//IzlLcoqFEPlJQBblxt1u46utAeyJdk3qoNcqhZZTvPfaVQxunRuM//jjHmbOXEVLTjHzRz31rq1N+IJtjmDszN2XirhkA9rs9V+TcobL+g8wao9CVu77JKdYCOFMdlmLclOaNYCduds9nfe5Sx7qiLdOg5+XFm+dhlG3NXfser6vZQQPX7/Ccf2BA3cyY8Ya2nv9y4ZD1/Dofa0do3t33H2pOBx/hTRTBsm65Zz3nmgPxtly6hT/8PAPEoyFEA4yQhblqjRqADvzdJOYu+dOvLMNJ8+8iyHxv47r/vqrF9OmraWTTww/HwshsIk93aigXGJvXf4vFUnpRqZs/IyL3h9h1VxytGsVHZNvl5xiIYR7EpBFuStpDeC8irtJLO9zs64sx5AY7ngdFXUHr722ji6+Z9h07FrqNPZ1eW9OFSi9RoPJauXJbtdy23X1adckt9RiTEoMo9eM45x2k8uz/7+9Ow9vusoaOP69WbpTlkLZF6UgUqSIjMhUGcRxF3VGZ8YR33EbGRVUFCgqIiCjA5TFBUQRRVTUURkFxI2xboCgoIWhiIAouywFCt2SJrnvH2nSpE3apFuWns/zvK82+SW5Gduc3HvPOTdB9+GVPyzgur6D6vy+hRDRSQKyiFh1SRI7cOAFduy42/3zli2DeOihlQxI2M0HP6bRrHXVx3vOstfvymfmJz/yyte7sWvNE384k50lb1U5p9igk2lZdhsphosZcnr/enjXQohoJQFZRCxfy8gWu4PEmKrdsDwdPLiI7dtHuH/eunUg48d/xLmJu3n/x55YYjSb9p7wuaSekhTLR1t+5YkPnfXDVrudUkMef11+F1a12+va5o7LaM9taENSveyVCyGimzQGERFtee5+spZuBqC0zEGsUaEMyu9e8q+/vsq2bTfjKoXatm0AY8euIk3/xIQPWxOTbKp2Tzq/0MKgaTlYbQ7sFHDc/DJFplVer+GqKT6j5TkN2pFMCBH+pDGIaDKu7teR3u2TueLprwCw2DXYNePe2USLBDPpHZq7g+GhQ6+zbdstuILx9u1nk5X1MWnspGDEIR755BCgsNj870nvO16C2ajJ159wwvwyDnXSPRZfNcUSiIUQgZKALCJekdVOrMmI1V5xXKHFprnzte9wlPeoPq/DWn744W9QHmh/+qkv48atorveyfE7DkOswqgMzp6aHgxKkXeggME9UwE4afuJn9U4SmK2el03IPVS/nPjC1LGJISoNQnIIuL5K0kqLj9v+PUv55Pcdxrg/Pnnn9MZM+a/dNc/kT/iCMrsjMJ27QCtqjzH3xd/yz//0JMNx19wnlOsKgK/Sady3zn/Yuaw2xro3QkhmgppDCIinmfDkQSzd0LX2anruKNPRTDevbsXY8Z8ygUtfmXCJ+2JT/A+pjH7+r7Emrz/LE7odQx/P5PstdnYHM5gbDKYuK3vaHbdt02CsRCiXsgMWYRUXY5i9HysqyQp78BJ7nhlAxabg4w23zKy3zSMBmcw3rOnJw88kMPgFod464eziIlVDOmTWuX1DUox6o3vsanDHDM/T4lxvdfrDu46mGeveJb01HR3lzBJ3BJC1JUEZBEydTmKsfJjJ17Vmz4dmpPeIZns6/vy8mcvcXfG45gMzhnt/v3dGTMmhyGtjvBm3lmYY5xL076alCTEagpM71BgegOtKtpxGnUyo/pPYc6w+1BKheQoSSFE9JKyJxES+YUWMqfnUFpWsfcbZzawZvzQGmeavh4LkBhjxK41s64uIKn4FrR2BtMDB05j9OgvuDDlGK9t6usOxr58tfsrRqy4k2353klbSbZLaVF2M4nmFqwZPxSg1uMXQjQdUvYkwl5dumz5eiw4s617tdqM+eRktNEKwK+/duWBBz7j96nHeeX7vpjMvoPx0eKjPs8pNju60apsJHGOM73G6Pr3xj9KUggRrSQgi5Coy1GM/rKqe7bcwv39pxBTHowPHerMAw/kcEm7EyzamMGJEgt5P58EtLs+2aEdvJz7MlmrssgvyXc/V6I5kXsGPMw7X55FmaMiyctzjKE4SlIIEb1qDMhKqZeAq4DDWus+5be1Av4NdAN+Af6stT7ecMMU0cbXYQ0jh6QF/VijQVFksZPW4gfGnDOJWJNzmfrIkQ488EAOl7YvZOE3Gbz/v/2MeSsXW3kMNRsV914ay1s/TWb1ntVez+95TvGg1P3uMbr2iV0zYM/xV75PCCGCVeMeslJqMFAIvOIRkGcAx7TW05RSDwIttdbjq3sekD1kUVV+oYUl6/cw77OdxBiDS45yZVm/v2EF/ZP+QbzZuZScn9+O0aO/4NIOFhZ8fRa7jpzi8qe/oszu/F13UEqB6Q1Omt4DZXc/X9fmXZl7xVyu6nmVz9fxlUldlyxxIUT0q9c9ZK31l0qpbpVuvgYYUv7vi4HPgRoDshC+PPv5Tiw2BxZbzUcoekpJiqW4aCP9k+50B+Njx1K5//7PuKRjCQvWZrBi837Gvb3JHYyLDes5ZvY+p9hkMDF2kP9ziqs7LrK+jpIUQoja7iG31VofBNBaH1RKpfq7UCk1AhgB0KVLl1q+nIhWtUnucs1KU2K28dMPlxNvdh53eOJEa8aMyaGtOsyof/fleLHzvGSrXZfXFC+gxLjO67kGdTqfF4Y9R3pqesO9SSGECECDJ3VprRcAC8C5ZN3QryciS7DJXa7a327NfuGeflkkxhQCUFDQijFjPiXFcYTDNxTRLM7IvuMlmAwOnzXFBp3MqP6TeXLYaJTyXwYlhBCNpbYB+ZBSqn357Lg9cLg+ByWajsrJXb6So1wz4sQYI+OXbiYldhej+j7kDsanTrVg7Nj/0tJxjCM3FBFnUhRZ7ewr+o4daiRWs/c5xZd2Hc4zV86gR5sOjfpehRCiOrUNyMuBm4Fp5f9cVm8jEk2Oq+2lZ3KUKwhv2V/A1JVbMRsMWGx2OiTtZczZD5MUdwqAwsLmjB27imTHCY78uRilwG44SfY3o1myZbHX6U0xuhuPDX6S8UOvCdE7FUII/wIpe3oDZwJXa6XUPmASzkD8llLqdmAP8KeGHKSIDtVlJHsmRy3L3U/WO5swKEVJeSesUhy0S9zHA/3G0yzOeQZxUVEzsrI+5tx28M0FpTQzGsjXn3AsfjFLtlRU4SWaE/lHvwcZlzmads2TGundCiFEcALJsv6rn7suquexiAgTTMlPoH2f8wstHvXCFSkHqQkHeOicB0mOLwCguDiJ8eM/4jftjLz86QDW7P6Oez8aya+H1oG14vk8a4qFECKcSacuUSvBHKyQX+jMdi4tc7izqV2lTYBXUM87cNLdvMOlTfyvPHzOgzRPOAFASUkCDz74AYnKwteZh7j+9ZdYvusF99GIAJ2Tu/DslfOq1BQLIUS4koAsglZdgPU1U/ZX2rTwq128uOYXYowKm0Mz47q+tEgwez22dfwhJpyTRYvEYwCUlsbz8MPvY3JY+fnqLzhmeI4dOytqitFGWnMdccduwFF8dgO8eyGEaBiGmi8RwpsrwHryPHShMl+lTcVWG/O/2IXV5qDQYqe0zMHYdzZT5jE9bhV3hIf7j6dFkjMYWyxxTJiwAm0oYdsfH+FI7FSvBh/xug/tLU+TWPo3rGUxZC3dTH6hBSGEiAQSkEXQgq0ddpU2xZkNNIs1EWtSXtnPLlabgzte3QhAi9ijPNI/i1bNjjrvs8YwceJ7tO5+lM3D/ujV4MOok3lw4NOcrrOJ0V3dt1f3JUEIIcKNBGQRtMoBNs5sqPFghav7dWTN+KG89veBvPC3AcSZjD6vc2hoHnuMif2zaJXsnP2WlZl59NF32dv+Gz5Jv8GrwUdzx2W8fMUaxp4/ApvDu++MnL4khIgksocsasWzdjgxxkiR1U5+oaXaoOwqbcovtGD3c6hJcsxxJvbLIqW5s9eMzWZi0qSlrG/1IfSf676ud+s+ZA2cyVW9hsjpS0KIqCABWdRaSlIsq3ceDTjb2vNxruBpAIrLa42bmQuY2G8crVv+CoDdbmTKlLf4uuWHcO58wFlTPGXIFO4deC9mo3cCmK8GI0IIESkkIIsqAq0vDqacqTLP4LnlQAGzP/6ah9LH0bpVRTCeOvUNVrf4EAa8ADhriidfMANtb83JEgcpPnp8yOlLQohIJQFZeAmmvthfOdOS9Xt49vOdNT6HK3j2bueg28lHKTMfAMBuN/DEE4v5ovmH0H+R+5xie/HZ/Hn+ZsyGn4I6N1kIISKBJHUJN88Z7ymLjdIyR7WlQ76yrS02B3NzdgT8HGVlJ/jms0GUmXcA4HAopk1bRE7yJ3D2K7Ry/IkZ569iUIeLgxqbEEJEGgnIwi3Y+mLPbOs4s/NxWmusdu+ELX/PYbMV8OlH/SmL2e6+LXvmC/w3+WNiz9pIe8vTNLPczKPLdpJ3oCCosQkhRKSRJWvh1qllPKU2u9dtpTZ7taVDV/frSO/2yVzx9FcAlDmqZk97lh+59qdbJ5awfs05tGp20H3dzFnzWZXwJR3S22Gy/BVVXqzsDMQqqNpnIYSINDJDFl50pXKkyj/7UmS1E+ujrjghxuhVo7wsdz+Z03O4YdGzLPu0h1cwnvPkM6wx7Kd7+pUk2i5yB2NwBt70DslB1z4LIUQkkRmycNt3vIR4s4lTlopDGuLNJvYdL6k28HVqGU+x1eZ1m1HBjOvOIjneTHqH5s5TnJZ+xXHDCzxy5hp6t6nY+336mTnkWpvTvO9pKIeBR4f1dp+B7FlPLGVNQohoJgFZuAXbEtPloy2/UmnbGA2MeXszMUYDVoeNc3pt5oD5Cab2snG2RzCe/+x08op6o/qVAXB1RgeGn9eVy/q08xl4paxJCBGtZMlauNWmJWZ+oYVHl22pcrtDOzOu86072a2yeOenh5l8hp1zUkvd17zw/ONsLjoHS3kwBnhrwz6WrNtNSlIsGZ1bSPAVQjQZMkMWXoJdFs47cLLK7BjAQSkFpjc4aXoPs8HOYz0TOLdtsfv+RQun8O2pgVjPKq3y2Ckr8risTzsJxkKIJkUCsqgiuGXhqtG42LCeY+bnsBuOYFIwKS2R89oVue9/ddEjrD2ZSVmfqsHYpaZ9ayGEiDYSkEWd7D1WUQdsU4c5Zn6eEuN6wJnYNbF7EpkdCt3XvP7KeL4oGIKttzMYmwxg8962xmrXJMb4Pg1KCCGilewhixrlF1rYtPdEla5Y+YUWpq7cisZGgWkpB2LvcgdjA/BotxYM7lgRjN96fSynnzUC01lW9x71fRf1JNbk/WsYZzZQZLVX+9pCCBFtZIYsqlVdb+t9x0uwGrZyMPYZygy73Y8xAE937Ux6l73u2/7z1r385ZqJ/O7KZO4o7OjeowaY9/nOKq/bqWV8UH21hRAi0skMWfhVXW/ro8VHyf5mND8bxnoF4zjdlXe7ZZDerSIYL1t6N5dfMokWfR3uM5NdGdT+MrsB6V0thGhSZIYs/PJ1mpPJAHPXv8AzGyeTX5Lvvl3pOFo7hrO4w2biu37rvn3leyPo2Gcc929cj/E7RZndwaRh6V51xr4yuzftPeHzJClJ9hJCRCsJyE1IoOccu3RqGY/VXhEQreoXdqn5TF6d53VdkmMQ7e0jeCTxKeJ75bpv/2jFbVxy6TTuWbuO0rKK55nw3hYeXbaFhBiT11K055hq26RECCEilQTkJqI2+7Grdx7F7nB41RSjKg6f6JzcBfuxWzFbz+HhpEl0GVQRjFd9+Deuu+ZpYk4rw/i1qvLcdo27RWfW0s1kprWu0pFrxnV9yao0ZpkdCyGilQTkJsBzL9i1BOwrCHraeegU497exEnWcyzWWVPsYjKYGDtoLCmOvzL3091Mjp9Ctwu+d9+f88mN/PHq+ZyTmUB+oYUyu8PXS7hVXop2zeQz01qzZvxQ6V0thGgSJCA3Ab72gqvbj12ybjcTVnzGYeNz7jIml/5tB/HKH1+gXUIav/3Xf5lsfoxuQza67//8v39i6CXzOCczAXDOdCcNS2fCe1Xba7p4LkVLZrUQoqmSgNwEBLMfu3jtTkZ/8DgF5jfQqiKj2aCTSXXczsfDn6B1szhyf8nnYcdUul66wX3N6s+vpazjYwwZ0sLrOYef1xUUTFmxFbNRYbU50FoTbzZ5LUXXZiYvhBDRQgJyExDofuzKbTmM+Pg2rObdXrcn2S4h1XErc/50Aa2bxaHtdo6+/me6Xloxe1775TBePXA3a8ee5nMMwwd25bL0dl71x5WXooOdyQshRDSRgNxEVHdoxNHio4xfNZ6Xcl/yqkw3O7rRqmwkzVRvPrz3AtLaNiO/oJhv5gwjfkiO+7pv1lzO4oOjyL67NwCb9p7wuefrGXg7tYwno7P3TFoyq4UQTZkE5Cak8qERDu1gce5ixq0aV6WmuIXtRprZrkZhYvIf0klr24zlG3aTv/gmTrtutfva77+5hF6DXmLtuS1ZvfMomdNz/O7/1rQ/LJnVQoimTGnt4+y8BjJgwAC9YcOGmi8UDS7vcB53rryT1XtWe91+XvvLOLz3r8SpVMrsDsZecgYdWsSz5/AJYv9zF93/8qX72u+/HUrfQW9w0fmp5BdayJye41VvHGc2sGb8UPf+cHX3ewq2XloIIcKVUmqj1npAINfKDLmJKbIW8dgXjzF73WxsDpv7dpNO5cHzZjD1spvdAXHL/gIeXZ6HKivjwRNz6D68Ihhv+m4wT/88BvOJDcxMyqBrSmK1+7/B7A8Hd/yjEEJEB+ll3YQs/3E5vZ/tzYy1MyqCsTaSXHY97Uuf5e01qe5e051axvPY+1sxWq2Mz3+atOFfuJ9nS24mT/2UhU7WWO2arKWbSYwxVrv/K/vDQghRPQnITcCegj1c++a1XPPmNewp2OO+PUH3ob3laVrabsFAnHvGCs7Eq1iblbGH59LjbxUJXFs3n8ecneNxNKt4frPBeVyir0MiXDNdf4dIyExYCCGcZMk6ipXZy5izbg5TvphCcVmx+/aU+BQmD/4Xz6zsiMUjh8BrRhsH9+x9lh53fOq+f1veb5i142HsSd7f41yPy+jcwm8mN1Sf6S2EEE2dBOQotXrPau5aeRdbDnt3yPpDz78x/ff/okebDnSJ218loxlg8/YD7JxxNz3u+MT9uO0/9GfWjxOwJ3oH4xijqjITri7Qyv6wEEL4JgE5zNQ1w9irpthD1+ReOI7dzo5tZzFs6yZmXKerzFhX7zzKxVNX8o+dL9Jj5Ifux+7a0Y/jrReimh2AirMliDEZ+OCe80lr2wwhhBB1IwE5jNSmj7MrgHdoEcuKnW+QtSrLq6Y40ZxI1m8nsuTT3ljKDJyi6glLrrKkx15fx99/XEyP+z7EYHAuZe/Z1ZdLr8ihY8eWnNkjtcqMWoKxEELUDwnIYaI2fZxdAdxu2MN+9Qwlyvuc4mt7XctTlz3F8ZPJvP3ZeixUlDlVLjk68Muv3Lb1Nc64f6U7GP+yqzddM96jY8eWgOwBCyFEQ5KAHCaC7eOcX2hh3NJvOKSXcFJ5n1PctXlXnrn8GYadMQyARKOl+pKjY8fInTieXvevwGh0XrfnlzOYvnUyK65s7fU42QMWQoiGIWVPYSLYOt3XNi3lF+OdnDQvrQjG2shtfUeTd3eeOxhDDSVHR4/y6q330n7kUoxG5/Ps29OD6XmPoeOTKLLafb28EEKIeiYz5DARaB/nPQV7uPfDe1n24zJQFbfH2vvQXo9kxsW3khhTdQbrc7n58GEW3zqWDqPfwmRyLmfv33c60/KmUGRMJA6kcYcQQjQSCchhpLo92jJ7GU+ue5LJX0z2qik26mTa6b+TqC8i+7qMwEuOfv2VF29+iM5j/o3ZXAbA4YOnM2fb46i4lsTJwQ5CCNGoJCCHGV97tP5qiv9+9t/JGjSF4tKE4JKsDhxgwf89SrfxbxATYwUg/3BXOvRZzptDu1BktUvSlhBCNDIJyGHMX03xWalnMf/K+WR2yQz+SffuZf5N/+T0h5cQE2MB4NjRTmRvmUpR3j7KHHuYcV1fOrWM93uusRBCiPonATkMObSDl3Nf9llTPHnIZO4beB9mozn4J969m2eGTydtwmvExpYCcOJYR7K3PM6v1lZQXhY15u1NGBTEGI0B10MLIYSoGwnIYWbL4S3ctfKuKucUu2qKuzTvUrsn3rWLOTc+Sa9JrxAf79yDPnmiPalnrKBw02HwqFEuszvrkC22qk1EhBBCNAwJyGHC3znFXZt3Ze4Vc7mq51U1Poe/tpt6+w5mDZ9H7ymLiI8vAuBUQVsyL/gcYrtS5vi12uetrh5aCCFE/ZCAHAaW/7icez68x+toRJPBxNhBY3lk8CMkxiTW+ByebTetdjujLuzBjQO70GrfL0y/cQF9HltEQkIhAIUn25B5wee0atUTwKvcqthqw669n7su5xbXtTe3EEI0FRKQQ8irptjD4K6DefaKZ0lPTQ/oeXy13Zy1ajsf/DuH877+if5PvERS0kkACgpa4Uh9k1aterkf7yq3yjtwkjte2YDd5t2g5NbfdqvV+6tNb24hhGiqpFNXCJTZy8hek82Z8870CsatE1rz8jUv8/nNnwccjKGi7aanMw79TL8vf+HsxxfSrNkJAE6dakF27lQe+sBCfqHF6/qUpFiax5uJMVb9lVi8djeZ03NYnrs/4DF5fkk4ZbFRWuYga+nmKq8rhBDCSWbIjWzNnjXcufJOnzXF034/jZSElKCfs3LbzTN/3cVZaw5x3swXSE4+DkBhYXNm5D7GvpKuNIv1vSfsq30nQHGZs31mMMldwfbmFkKIpk5myI3kaPFRbl92O+cvOt8rGPdJ7cPqW1fzwtUv1CoYQ0XbzViTos+BHfT+Mp9BMxbQvLmzZKqouBnZ309hb9HpgP89Yc+e1wkxxir3uwJqIILtzS2EEE2dBOQG5tAOXvr+JXrN7eXV4CPBnMCM38/guxHfVWnwkV9oYdPeE0Et717dryPrhiRxxupCMmc+R8uWRwAoLW2GrdXrHLL2rHqwhI/XzExrzZrxQ3nupv7EmpTXNcEE1GoPtBBCCFGFLFk3oNrUFNc2Ecrx5Wqm/n0tv535LK1aHQLAaknkjLNWcNKRwfujjH5bYvp7zezrM2o87KI6cn6yEEIETmmta76qngwYMEBv2LCh0V4vVKqrKfY8p7iy/EILmdNzKC2rWOqNMxtYM35otcHMLp7mwwAAIABJREFUnvMF99/yDRfMeoo2bZyJV2XWREqav0rW+3HVBveaXlPKloQQovaUUhu11gMCuVZmyPWsNjXFrqBXUGKtkghlVIrPth3mwl6pPgOi/ZNPuefW/zF41lx3MLaVxZPW+10uf97mVQrlKymruuQr1/0SjIUQouFJQK4n/mqKL+hyAfOvnO+3jKlyQw9HpQWLIqudySvyeGTZliozXNsHn3DX7TsYOvsp2rVzfgGw2eI4u/8K9hWfg9mwvsYsZ3/JV1v2F/CXBV9LDbEQQjQSSeqqo+pqihdds4gvbvnCbzCuXKtrsWm01sSaDCR6ZDkXWuxV6njLln/IiFt3MWTmbNq3/wUAuy2Wfv2WkZJyUcBZzr6SryZe2ZupK7cGXUNcm2Q0IYQQTjJDroO61hT7Wi6ON5uYN7w/h0+WMnlFHoUWu/s+1ww3+eOP+fudh/n97Nl07LgLALs9hox+/0HF/c59bKJnS8zqkrIqJ1/VpoZYunIJIUTdSECuJJAkpvzifMb/dzwvfv+i1+3BnlPcqWU8VnvVWWx6h2TSOyTzyLItVe7rlPMJtzxQxKVzZtK58w4AHHYzffu+w5p9GYxfmuMVFNeMHxrQPnBKUqzX/cHUEPtq3SknRAkhRHBkydrDstz9ZE7P4aaF6322inRoB4u+X8QZc8/wCsaJ5kRmXjyTjSM2BhyMAVbvPIrdI/CZjco9i61o9uFs0hFrMvBS3B7uHm3h4lkz6dLlR+eYHCbO6vtvDPGX+GxVCZDRuYXfumNfy8vB1hD7at0ZTBMRIYQQMkN2q2mWVx/nFHvOvgHGL92M5zkOBgWZaa3dP2vX/9eKK3K/YPZXGVw++19067YVAIfDSJ8+b9CmzR/YtPdEwMvMgSwvB1NDLF25hBCi7iQgl/O3b7rj8FFmfP100DXFlVUOgiOHpFV5vRij0R1AXV8QLDbNtd9/yq6vz2XY7Mfp3v1/AGhtID39NVJTrwcCD4rBLC9XXsb2xzWjrksTESGEaOokIJfzFdCO67X86b072Hdqr/u2YM8pBt9BcO5nOwD/rSn3HS/BqBTXbfyU7evOZ9isx0lL2wSAw6EoTXyStm1vcD/WV1CceGVv97KxKzg21KEP0pVLCCHqRgJyOc+AZleHOajmU2RYD6cqrqmpptgfX0EwxmhkxODTmff5Tp+zyi37C7hyzUf8b/1Qrs7+Jz17fgc4g/GLW0az8UgaA9MtXoHPMyhu2V/A1JVby+ubHYy6MI0bB3Zp0OXlQGfUQgghqpKA7OHys1J544ev+PePs9GqItmpdUJrsi/O5uaMm1FKVfMMvvkLgjcO7MKNA7tUmVXmF1r4ccIsctddxjUzpnLmmd+6H7co7x7WHLjI7xGKrp//suBrrxn5rFXbmfvZDrKvz5DlZSGECEMSkMut2bOGO1b8gx+O5nmtJDd3XMraW16iR5sOtX7u6paTO7WMJ6NzC6/r86c8yZq1V3LtjKmkp69z3/5y3ki+2n8JUP2s1teMHMBi02Qt3cya8UMDLocSQgjROOoUkJVSv+Bc1LUDtkAbaIcTfzXFZkc3WpXdTRtzX4pLE+r8Ov6Wkz2znPMLLeQ/MpM7nj+PP057jLPOWuN+/JIf7uLrg1fQLNZY46zW14zc/b7K94t9lUIJIYQInfqYIV+otT5aD8/TqBzaweLcxYxbNY78knz37UrH0cJ2I81sV6Mw1Wv5jr/l5KylmzlVUsbesf9i5ca/8ocnppKR8aX7cW9vv4NhgyYwLcCkKdeMfNw7m7HYpBxJCCEiQZNcss47nMedK+/0WVN8VeeHyP4wH3Nsw+yv+lpONio4OHoqy76/lT8+/hj9++e474tPeYzHbsxyjyHQsbhm5K+v38Pcz3YQY6x5Zi2EECJ06hqQNfCJUkoDz2utF9TDmBpMoOcUX9u34c4ArrKcrDW3f/gWb2z6B9c99hi/+c0q911xrSYy8KyJtX6tlKRY7rmoh8/EMSGEEOGlrgE5U2t9QCmVCqxSSm3TWn/peYFSagQwAqBLl5q7WTWUFT+uYNSHo6qcUzxm0BgmDp7oVVPsWb4TSG/ryqp7jFeCl1KMWPk6L2++h+sfm8J5533ovm75TzcxYcCEurxln+9HCCFEeFJa65qvCuSJlJoMFGqtZ/q7ZsCAAXrDhg318nqBqu05xVC7E4wCfUz+qVKOjHiYG94ezvWTp3L++RXjW/nzX7ng7Gw5LUkIISKcUmpjoAnPtZ4hK6USAYPW+lT5v18CPFbb56tvZfYynlz3JJO/mExxWbH79pT4FLIvzuaWfrdUW1NcmxOMqnsMUDFrTjBjGzmRP7/9f1w38XGvYBzX4n7GX/8ErZvF1fl/AyGEEJGjLkvWbYF3y4OaCXhda/1RvYyqjup6TjFUtK70ZFCKvAMFpHdo7nNJ2l9byiXr9/BseUcum93GKxvf55ZP7ua6CY/zu98tdV/bufM4Tj99ut8vCrVZPhdCCBEZah2Qtda7gIx6HEud1dc5xfmFFj7acpAiq93r9mKrnVsXfYvBoIgzGassSfuq/7Xa7cz7bCcWmwOro4yHlr/G/+2YzHUPPsHQof92X9ep0+hqg3Ftls+FEEJEjqg4D1lrXW/nFC/L3c9vp33K/C92+bzfrqHMrr3OHHadKezrHOFRF/YgxmjA6LAz4d1Xyd4+lWvGZnPxxUvcz9mx4yi6d59d7cy48lnH497Z7PMsYyGEEJEp4uuQ8w7ncdfKu/hqz1det19zxjU8ffnTAZ1T7OJ55GGgjAbl1VO68qlHAM9/uo0JSxfzxM/TuPaBGVx++cvux7dv/w/S0p6udj/bdWKTJ4vNwevr93DPRT0CHqsQQojwFbEz5CJrEeNXjaff8/28gnHX5l1ZfsNy3rvhvaCCMVTsAQc1DoudLfsLAGdA37T3BIC7NWVKjOKlz99h6q5sht03m6uuWuh+bLt2t9Oz57M1HliRGGOktKxqK8xncnbILFkIIaJERM6QV25fyd0f3B1QTXEw/PWANhsVWmvizSaKrTbslSbQU1dudf/Ta4/3zNbsunwkf/p6CtfcM4trrnnO/Zi2bW/mjDMWoFTNXwCKrHbMBkWZw/uFTca6n2MshBAiPERkQN56ZKtXMK7tOcWVVT6VyWq3M+pCZ6crcM6gC0qs3PXad14JX0almPL+Vqy2inKniW9+S9rqZVy6ZgpX3z2HP/5xrvv61NTh9Or1IkoZAsqc7tQyHoMBKh3ehN2hpS+1EEJEiYgMyKPPG80rm1/h18Jf63ROsS+V94A9g2RKUiz5hRbslZqplNkdxJgMWMu7ccaVlfLwW69y8b65DPvHU/zpT0+6r01NvYFevV5GKWPAmdMpSbFkX5/BA2/l4jorwmxUZF8vfamFECJa1FunrkDUZ6euH478QGpiakA1xbVR3cx1ee7+KmcbT1q+BZsDEqwlTHzzFcYcnM+Vf5/L8OHT3I/beOh8evR8javP7kp+oYXM6Tlee8NxZgNrxg+ttvFI3oGTgCa9Q3MJxkIIEeYapVNXqJ3Z5sygrg+mqUZNM9fMtNbMvD6Do4UWzk9rTcvEGCavyCPRUsTEN15l9KEFXHHrs17BeMOvg5i/aSzmLVvJ7NHObxOR6vaEU5JiGdyzTVDvWwghRGSI2IAcjECWhl0BOzHGWG3LzGW5+xn79ibKyjO7TAa476KeJJcWMX7JK4w68iKX/+05/va3qe7n/v7wuczflIVdm0goD7q+EsgCOatYunUJIUR0ivqAHEhPas+AbbHZMRi896ONSvHZtsP069yCrHc2u4MxgM0BL7//HVmvLWHk0Ze45MaF3HrrJPf9W44OYN73D2HXZqAi6FZOIAvkrGLp1iWEENEr6gNyTUvDvgJ25bqmIqudySvysNocgHewblFykodee4N/HFvExX95iTvueNh9nyl+KF26L8C0aTvxPoJudQlkldXmsAshhBCRI+oDck1Lw74CdqxRoZXCbFDu8qZCi6vMqSJYtyou4KFX3+SOE4u56LrF3Hlnlvu+H/IzuP7ipbRJbkFmj07Vno8cSECtzZ4zyBK3EEJEiojt1BUoX/2lPWepvgK2Mig+uOd8plydTlKs0eu+WKPCZFC0LjrOQ6+8yR0nXmHItUsYNep+9zXbjvXhmdyJGAzx7jG4One5uLp6BdppqzZ7zsty95M5PYebFq4nc3oOy3P3B/RaQgghGl/Uz5Ch5tpiX3u5aW2b0TIxhkeWeR/hqAyKVdd3Z/cVE7im4FV+N+wN7rvvHvf9Px7rzZyNkzAZE/3OXmuzFxzsnrMscQshRGRpEgEZql8a9hewfXXuykpPYu8VE7n6wItccMVbPPDAXe7n2XG8F3M2TsZij0fje/Zal0AZzJ5zbZe4hRBChEaTCcg18RewXUFwyfo9vPveGpLvWMAVha9x/qXvMGbMCPd1u070YPbGKZTaEwCYNCzd5/MFGij97f0Guudc27IqIYQQoSEBOUDL3l3N3YuWc3PxEn77+3fJyroNg8GZ4PVzQRozN0ylxOY81MKo4LI+7Xw+TyCBsj7Km2pTViWEECJ0JCDXIL/QwjervuEfi97n5uIlnHfhch588GZ3MC7RZzJzw6MU25Lcj0mIMfldGq4pUNbn3m8wS9xCCCFCK2oDcn2U+yzL3c9zCz7g1kUfcXPpawz83ftMmHATRqMzUMbF9+HMnh9h+3Qznkcx1bQ0XF2grO+930CXuIUQQoRWVAbk+ljyzS+0sOC597n55f/yN8sSBmR+xCOP3IjR6KxHtht60v/sT4mJSSX7eoJeGvYXKGXvVwghmqaoC8j1seSbX2hh4/tfMvzlHG62vMY5gz5h0qQ/YzKVn69oSiM9/RNiYlKB+l0alr1fIYRomqIuINd1yXdZ7n4WzXuXP7+ymlusr3L2uZ8yefL1mM1lABwq6sDTuVPI/3grM64zuWfe9bk0LHu/QgjR9ERdp666LPnmF1p49Zl3uH7xWm62vkrGgM+ZOvUPxMRYAThc3I5p3zzB/lPNKS1zkLV0c8Cdtvy9nr9uXb66ewkhhIheUReQa2qVWZ1jOV8x7JUN3FL2Cn37f8k//3kNMTHlwdLYmbm5Mzhuae2+3jXz9hRoS0xpaymEEMJT1C1ZQy2XfNes4du/LOF22yL6ZKzh8ceHERtbCoA5pjOnn7GKw5/sorps6kCTyaStpRBCiMqibobsEtSS7xdfsOSiF/m/0oX0Putr/vWvK4mLc858Hao9/c/+jPYpZ1Q78/YMsqcstmqXtF373J58zbaFEEI0HVE5Qw5GwcqP+M8flnJ72UJ6p69j2rQriI8vAsBkbk//s78gPr47UH/1w1LaJIQQorKonSH747nHu3b+67x59X+4vex5zui1genTLyMhoRAAs7ktZ/fLISGhh9fj/R2lWFBixWq3e13rL8jWZZ9bCCFEdGpSM+Ql63YzZUUeZqOB325fT/+3DnK3XkDPnhvIzr6ExMRTAJjNbejXL4fExF41dvzy3Dd2aDAZIN5sqrF+WEqbhBBCeIr4gBxoi8wl63Yz4T3n2caDt60l49187tbPkZb2PdnZl5CUVACAMqSQkZFDYmLvGpO0fCVnxZoMzBven/QOybXu1iWEEKLpieiAHExW85T3twJw2Y9rOHNZIffoeZx++mZmzryY5OTjABSWNaNfxockJfXxGWzHvbOZFgkx7mDra984xmigebxZAq0QQoigRGxADqZ0aN/xEmKMikv+9yWnLbcymqfp1i2PWbMuonnzfACKypKg5Vt0Sv2N+zGVg63F5uDOVzfiQDPjur5kprWW5CwhhBD1ImKTuoIpHerUMp7LNufQdbmdscyhS5cfmDVrKC1aHAWgxJbAU99PxWE6C3AlaZVhtTuqPFdxmd1d0gRIcpYQQoh6EbEz5GBKh1KWvknr9wxkMYNOnbYza9ZFtGp1GIASWzwzv32Mnwp6kLV0MwcLSpm1ajsxRoXd4cBsVJiNBoqt3hnUruAvyVlCCCHqQ8TOkAMuHVq4kMdv2cFDTKNDh5+YPXsorVsfBKDUFsesDVP4qaAXAFpr/vXhNqw2B4UWOzYHKGDkkO7EGJXX03oGf+k7LYQQoq4idoYMAZQOzZ/PlLt/ZTL/pF27n5k9+0LatHH2jLbYY5m9cTI7T/R2X26x6SqvYbVr5n32Ew6tMRsVcSajV0lToFneQghRF43xWSOfZ6EV0QEZ/JcO6aeeZtLoE0xlCm3b7mbOnAtp23YvAFZ7LHM2TmL78T7u680G0Bp8xGSKy5zL1bEmmDf8bNI7NCclKTbgLG8hhKiLxviskc+z0IvYJevq6JmzeGT0KabyKG3a7GXWrKG0a7fbeaeKZcH/JrHtWF/39WaDAqV8BmNPMUYjzeNj3DPjQHtXCyFEbTXGZ418noWHqAvIhZP/yYPjbDzBBFq33s/s2UPp2HEXAErF0C3tbbbk9/N6TJlDU2b3jsYmA1TaNvbaN5YDIoQQjaExPmvk8yw8RFVA3nrXOCZNiWcG42nV6iCzZg2lU6edAChlpk+f/9Ct4zCvZLAYoyLO7P0/Q0KMkZdu+Q1z/tLPb9KYHBAhhGgMjfFZI59n4SHi95AB0JqiByfw/HOdeJr7aNnyELNmXUSXLtsBUMpEevo7pKRcCXgngyXGGLlq7mqvp3No7d4n9pc05sryzqq05yKJEEKI+tQYnzXyeRYelNY1bJzWowEDBugNGzbU75NqjeOhCdwzvSPPMpLmzY8wZ86FnHZaXvkFRtLT/02bNtd5Pcwzm3DNzqNVfhEDTWaQrEQhRGOQLOvIpJTaqLUeEMi1kT1D1hrHmHHcPSeN57mT5OR8Zs36vTsYO7SBEzFPeQXj/EILS9bvYd5nO4kxVgTgNeOHev0iBvqLKQdECCEaSuXPoYb+rJHPs9CK3ICsNY57R/OPuX1YyB0kJR0nO/tiund3trR0aMWCzQ+Qe/R0fpdhcZcpZb2zyV1vbLFV9MBeM34oGZ1bAJL+L4QIPfkcanoiM6nL4cB+50hun9uPhdxBYuIJsrMvoWfP7513a8XC/41m3cEh7kxBV1q/r+YfntmEkv4vhAg1+RxqmiIyINvvH8utC87jZW4lIeEkM2ZcRq9eFXvTi7bcw9oDFwEVmYK+0vpdAi1nyi+0sGnvCfmjEEI0KClDapoicsn6vgNZvEo74uNPMX365fTuvd59X0n8NL493JdmsVUzBSun9QPEmgIrZ9qyv4C/LPhalo+EaOIaI/GpujIkSbyKXhGZZZ2bC5dffIIHH7majIyv3Lf36DGXjh1Hev3CAj6zqa12O6Mu7MGNA7tU+aVenrvfK+t64lW9mfr+VkrLKv5A4swG1owfKn8QQjQhjbmvW/lzaMZ1fdEg+8oRJpgs64gMyHZ7MWvXXoXd/pn7trS0J4lvcadXIPaVTR3oUYmeQX3f8RJuWrieUxab+/5msSZe+/tAdyKYECK65RdayJye06hfzCtPLhr79UXdRX3Z044do7yCcffuM9l49HrGP5+D2WCg1GbH4dC4umH6y6YG/8tPldP/pYuNEE2ba1+3lIrPAte+bkMFRM/PoU17TzT664vGFZFJXV27Pkpc3GkAnH76NBJajvLKSCyzVwRjT5WTIpbl7idzeg43LVxP5vQclufu9/l6AZ+9LISIWnVtL1nXpFBpbxn9InKGHB/fjX79viA/fzkdO47ky+2HMShV4+M8f3k9ywpc3zizlm4mM621z0CbmdaaBf83AKhoqymEaDrq0l6yPvaepb1l9IvIgAwQF9eZjh1HVmn24U/lbOpglp+kQF8IAd598APNcg72y399v76IHBEbkMH5i571ju9mH55MBsXKe84nrW0z955xYozR5/JPYoyRTXtPeLXQrK8/JiFE5Au2vWR97z1Le8voFdEBecn6Pe6ELZdYkwEDUOJxe7zZSJHVXj6b3ozRoLA7NH/5TSfe2rDPPfP984BOXDV3tddMuGtKoiRSCCFqTfZ+RaAiNiDnF1p46r/bfd7nwHvG7Jr5jn17E2Ue2V5vfLOXD++9gCKr3X0MY+WZ8Pujzpc/JiFErcnerwhUxAbkvAMnfWZS35bZjTPbJ1f55T9QUOIVjAHK7JoDBSUM7pnqt6SgyGqXPyYhRJ3I3q8IRMQGZPC9bzyoewqDe6ZW+eX/cvsRP8/jzM6ublkpo3ML+WMSQtSJ7P2KmkRkHTJAeofmmI3epU5moyK9Q3PA+cuf0bmF+w8gvUMypkrv1mRw3u66vrpa48rPJ4SIHOF2MEx14wm3sYrGE7Ez5JSkWGb9KYNxHkla2df7X0pOSYpl9p/7MfbtXBQGNA5m/qmf1/WyrCRE9Am3ssXqxhNuY3WRAy0aR8QGZAg+gGpAKUN5APe9OCDLSkJEj3ArW6xuPEBYjdUlXL8kRKOIXbJ2cS0lA9Uu8+w8dIpxb2/CYnNQbLVjscmB30JEu3A7V7i68YR6rL6Wyj2/QJyy2Cgtk8/NhhTRM2SXmr7BLcvdzwNv5WKvdByy1BMLEd0aowY4mOXcmsYTqhJLf5+hoThQoymL+BlyTd/gXN28KgdjAKtd6omFiGaBHgxT20SqZbn7+e20HP76wjp+O83/ATWBjKchDrEJ5H1V9xkqTU0aV8TPkGv6BrfveAlGg++DJ0ZdmCbf8oSIcjXlmtR2jzS/0FKl2dCYtzfVuOdb3Xg870uMcXYYzC+0eF0T6Iw80PdV3WdoRucWXn0YrHYHI4ek1fi/jaidiA/INX2D69QyHrujas1yjBFuHNilUcYohAitysmanj3ta5tIlXegwGezobwDBQzumRrUeCrft3rnUZ/BNNAgG0wyW02foa4vCUvW72HeZztY8OUu5n2+s9rXlozs2on4JetA6oezr+/rVbNsMlCl5EkI0TR4noN+xdNfVbk/8EQqf0e+1nwUbHX8LSHvPHQq4ASrYBLEAl0qf/bznVhsutrXDvSMeV/vWWqvo2CGDFWXgACvE5tc9+cdKAAU6R2SJRgLESWCmZH5mjlW7sEb6B6pq9mQ5/k2ns2GasvfEnKun/a+vhKsgt37rWlZP5DkrtqWmElZVYWoCMiA+z+4c1llJzFG7/+4KUmxNS4jCSEiS7Af5r4CS6xRoZUi1ui7V72/gO9qNjTunU0YUNi0g8nD+tT5y76/YNqvc4uAg2xtDrSobhk9kABfm4zscKsTD7WoCcjOoxU3uc9Gdh3L2JT/4woRzWrzYe4rsCiDYslt5/JLfjH9OrcgrW0z9301Bfyr+3XkVKmNKe9vJdZkZOrKrTSLM9VphucvmKa1bRZUkK3PzoOBBPjaZGRLWZW3OgVkpdRlwFOAEViotZ5WL6MKkusP0xWMPTXl/7hCRLPafJj7Cix/PqcTN730TZWgG0jAzy+0MHXlVqw2B1YbPq+pDX/BNNggW5+dB2t67drMysOtrCrUCWm1DshKKSMwD7gY2Ad8q5RarrXeWl+DC5SvP0wXqZkTIjrV9sO8cmmRr3PQXffXFPAD/VJQmw96f8E0lO19a3rt2nxhCJfjbcNhL7suM+RzgZ1a610ASqk3gWuARg/Ivv4wAWJNdS+sF0KEp0A+zKvb/01JivV7DrrrMTUF/ECuCYcP+sYU7BeGcDjUJ1z2susSkDsCez1+3gcMrHyRUmoEMAKgS5eGqfut/IdptdsZdWEPbhzYRYKxEFGsug/zQAJhdQE1kIBf0zXh8kEf7kJ9qE+47GXXJSD7KrarsomrtV4ALAAYMGBA1U3eehIO37KEEI3P14d5MIFw5JA05laqzAhmz7a6a8Llgz7Ue6PhLlz2susSkPcBnT1+7gQcqNtw6ibU37KEEOEhkEDoOYMGzYjBp/tcVQvkc8XfNeHwQd/UlsxrI1z2susSkL8FeiilTgP2AzcAN9bLqIQQog5qCoS+ZtDzPt9Z7+10G+uD3t8MWJbMAxcOq6y1Dshaa5tSahTwMc6yp5e01nn1NjIhhKilmgJhYy4lN/QHfXUz4HBZMo8UoV5lrVMdstb6A+CDehqLEELUG3+BML/QQkGJFavd7nV9Qy4lN9QHfU0z4HBYMheBi5pOXUKI8BPqZKLKgdBzNunQzt7T8WZTSOtf66KmGXC47I2KwEhAFkI0iHBLJvI1m4w1GZg3vH9QB86E+kuGp0BmwOGwNyoCIwFZCFHvwjGZyNdsMsZooHm8OeAxhduXjEBnwKHeGxWBkYAshKh34ZhMVNf91HD8kgEyA44mhpovEUKI4IRjMpFrNhlnNtAs1kScObjWuq4vGZ5cXzJCLSUplozOLSQYRziZIQsh6l24JhPVZTbp70tGYoyRTXtPyOxU1JnSusG6WVYxYMAAvWHDhkZ7PSFEaIVTAlR9WJ673/voxgGdeGvDPkwGhdWumTSsN8MHdg31MEUYUUpt1FoPCORamSELIRpMtCUT+Tu60WXCu1tAw/DzJCiL4MkeshBC+JFfaGHT3hPkF1rct7n2a4usdkyGqmfsTFmR53W9EIGSGbIQQvhQU4lTp5bxWO1Vt/zMxsbJJo+27QAhM2QhhKjCs8TplMVGaZmDrKWbq8yUJw3rXeWxdq0bPJt8We5+MqfncNPC9WROz2F57v4GfT3ROCQgCyFEJYGWOA0f2JXHr+1DjFGRGGMMupSqNgL5siAikyxZCyFEJcHUUQ8/ryuX9WnXaMvH4dh0RdQPmSELIUQlwTYRaczGHOHYdEXUD5khCyGED+HakjJcm66IupOALIQQfoRrHXW4flkQdSMBWQghIlC4flkQtSd7yEIIIUQYkIAshBBChAEJyEIIIUQYkIAshBBChAEJyEIIISKOr4M/Ip1kWQshhIgoNR38EalkhiyEECJiRHMvbwnIQgghIkagB39EIgnIQgghIkY09/KWgCyEECJiBHvwRySRpC4hhIhw+YWWJtXXOlp7eUtAFkKICBYMSEoFAAAE6ElEQVStGcc1icZe3rJkLYQQESqaM46bIgnIQggRoaI547gpkoAshBARKpozjpsiCchCCBGhojnjuCmSpC4hhIhg0Zpx3BRJQBZCiAgXjRnHTZEsWQshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQAKyEEIIEQYkIAshhBBhQGmtG+/FlDoC7K7Hp2wNHK3H54sU8r6bjqb4nkHed1MTze+7q9a6TSAXNmpArm9KqQ1a6wGhHkdjk/fddDTF9wzyvkM9jsbWVN93ZbJkLYQQQoQBCchCCCFEGIj0gLwg1AMIEXnfTUdTfM8g77upaarv20tE7yELIYQQ0SLSZ8hCCCFEVJCALIQQQoSBiA3ISqnLlFI/KqV2KqUeDPV4GoNSqrNS6jOl1A9KqTyl1H2hHlNjUUoZlVLfK6XeD/VYGotSqoVS6h2l1Lby/+aDQj2mxqCUur/893uLUuoNpVRcqMfUEJRSLymlDiultnjc1koptUoptaP8ny1DOcb65uc9Z5f/jm9WSr2rlGoRyjGGUkQGZKWUEZgHXA70Bv6qlOod2lE1ChswRmt9JnAeMLKJvG+A+4AfQj2IRvYU8JHWuheQQRN4/0qpjsC9wACtdR/ACNwQ2lE1mJeByyrd9iDwqda6B/Bp+c/R5GWqvudVQB+tdV9gO/BQYw8qXERkQAbOBXZqrXdpra3Am8A1IR5Tg9NaH9Raf1f+76dwfkB3DO2oGp5SqhNwJbAw1GNpLEqpZGAw8CKA1tqqtT4R2lE1GhMQr5QyAQnAgRCPp0Forb8EjlW6+Rpgcfm/LwaubdRBNTBf71lr/YnW2lb+4zqgU6MPLExEakDuCOz1+HkfTSAweVJKdQPOBtaHdiSN4kkgC3CEeiCN6HTgCLCofKl+oVIqMdSDamha6/3ATGAPcBAo0Fp/EtpRNaq2WuuD4PwCDqSGeDyN7Tbgw1APIlQiNSArH7c1mfotpVQSsBQYrbU+GerxNCSl1FXAYa31xlCPpZGZgP7AfK312UAR0bd8WUX5nuk1wGlAByBRKXVTaEclGoNSagLObbkloR5LqERqQN4HdPb4uRNRuqxVmVLKjDMYL9Fa/yfU42kEmcDVSqlfcG5NDFVKvRbaITWKfcA+rbVrBeQdnAE62v0e+FlrfURrXQb8B/htiMfUmA4ppdoDlP/zcIjH0yiUUjcDVwHDdRNujhGpAflboIdS6jSlVAzOpI/lIR5Tg1NKKZx7ij9orWeHejyNQWv9kNa6k9a6G87/zjla66ifMWmtfwX2KqXOKL/pImBrCIfUWPYA5ymlEsp/3y+iCSSzeVgO3Fz+7zcDy0I4lkahlLoMGA9crbUuDvV4QikiA3J5AsAo4GOcf6xvaa3zQjuqRpEJ/B/OWWJu+f9dEepBiQZzD7BEKbUZ6Ac8EeLxNLjyFYF3gO+A/+H8jIrKtopKqTeAr4EzlFL7lFK3A9OAi5VSO4CLy3+OGn7e81ygGbCq/DPtuZAOMoSkdaYQQggRBiJyhiyEEEJEGwnIQgghRBiQgCyEEEKEAQnIQgghRBiQgCyEEEKEAQnIQgghRBiQgCyEEEKEgf8HrF6l0UHXy+QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from sklearn.datasets import make_regression\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(8,8))\n",
    "np.random.seed(42)\n",
    "\n",
    "class robust_regression:\n",
    "    def __init__(self,loss_func=\"huber\",delta=1.0,valid_freq=100):\n",
    "        '''\n",
    "        Class initialization\n",
    "\n",
    "        Parameters:\n",
    "            loss_func:string(\"l1\",\"huber\",\"l2\"), optional, default \"huber\"\n",
    "                The loss function\n",
    "            delta: float,optional,default 1.0\n",
    "                The delta parameter of huber loss\n",
    "            valid_freq:int,optional,default 100\n",
    "                validate frequence\n",
    "        '''\n",
    "        self.loss_func=loss_func\n",
    "        self.delta=delta\n",
    "        self.valid_freq=valid_freq\n",
    "    \n",
    "    def fit(self,X,y,X_val=None,y_val=None,max_iters=2000,batch_size=32,learning_rate=0.01):\n",
    "        '''\n",
    "        Train a robust regression model\n",
    "        \n",
    "        Parameters:\n",
    "            X:array-like,shape(n_samples,n_features)\n",
    "                training data\n",
    "            y:array-like,shape(n_samples,)\n",
    "                traing target\n",
    "            X_val:array-like,shape(n_samples,n_features),optional, default None\n",
    "                validate data\n",
    "            y_val:array-like,shape(n_samples,n_features),optional, default None\n",
    "                validate target\n",
    "            max_iters:int, optional, default 2000\n",
    "                maximum iterations\n",
    "            batch_size:int,optional,default 32\n",
    "                batch size\n",
    "            learning_rate:float,optional,default 0.01\n",
    "                learning rate\n",
    "        '''\n",
    "        with tf.Graph().as_default():\n",
    "            shape=X.shape[1]\n",
    "            self.weight=tf.get_variable(\"weight\",shape=[shape],initializer=tf.random_normal_initializer())\n",
    "            self.bias=tf.get_variable(\"bias\",shape=[],initializer=tf.constant_initializer(0.0))\n",
    "\n",
    "            X=X.astype(np.float32)\n",
    "            y=y.astype(np.float32)\n",
    "            train_dataset=tf.data.Dataset.from_tensor_slices((X,y))\n",
    "            train_dataset=train_dataset.batch(batch_size).repeat()\n",
    "\n",
    "            val_dataset=train_dataset    # The val_dataset must be defined,regardless of training or validate\n",
    "                                         # because of the use of tf.cond\n",
    "            valid_flag=(X_val is not None and y_val is not None)\n",
    "            if valid_flag:\n",
    "                X_val=X_val.astype(np.float32)\n",
    "                y_val=y_val.astype(np.float32)\n",
    "                val_dataset=tf.data.Dataset.from_tensor_slices((X_val,y_val))\n",
    "                val_dataset=val_dataset.batch(batch_size).repeat()\n",
    "\n",
    "            training=tf.placeholder(dtype=tf.bool)\n",
    "            x,y_true=tf.cond(training,lambda:train_dataset.make_one_shot_iterator().get_next(),\n",
    "                            lambda: val_dataset.make_one_shot_iterator().get_next())\n",
    "\n",
    "            y_pred=tf.reduce_sum(tf.multiply(x,self.weight),axis=1)+self.bias\n",
    "            loss_l1=tf.losses.absolute_difference(y_true,y_pred,\n",
    "                                                  reduction=tf.losses.Reduction.SUM_BY_NONZERO_WEIGHTS)\n",
    "            loss_l2=tf.losses.mean_squared_error(y_true,y_pred)\n",
    "            loss_huber=tf.losses.huber_loss(y_true,y_pred,delta=self.delta)\n",
    "\n",
    "            loss=tf.case({tf.cast(self.loss_func==\"l1\",tf.bool):(lambda:loss_l1),\n",
    "                         tf.cast(self.loss_func==\"l2\",tf.bool):(lambda:loss_l2)},exclusive=True,\n",
    "                         default=lambda:loss_huber)\n",
    "\n",
    "            train_op=tf.train.GradientDescentOptimizer(learning_rate).minimize(loss)\n",
    "\n",
    "            with tf.Session() as sess:\n",
    "                sess.run(tf.global_variables_initializer())\n",
    "                for i in range(max_iters):\n",
    "                    loss_,_=sess.run([loss,train_op],feed_dict={training:True})\n",
    "                    if i%self.valid_freq==0 and valid_flag:\n",
    "                        loss_l2_=sess.run(loss_l2,feed_dict={training:False})\n",
    "                        '''if i==0:\n",
    "                            if self.loss_func==\"huber\":\n",
    "                                print(\"The loss function is %s and the delta parameter is %.4f\"\n",
    "                                      %(self.loss_func,self.delta))\n",
    "                            else:\n",
    "                                print(\"The loss function is %s\"%self.loss_func)\n",
    "                        print(\"The %d th iteration loss is %.4f \"%(i,loss_l2_))'''\n",
    "                weight_,bias_=sess.run([self.weight,self.bias],feed_dict={training:True})\n",
    "                self.coef_={\"weight\":weight_,\"bias\":bias_}\n",
    "        \n",
    "                    \n",
    "    def coef(self):\n",
    "        return self.coef_\n",
    "\n",
    "n_samples=500\n",
    "X=np.random.rand(n_samples)*10.0\n",
    "true_w=2.0\n",
    "true_b=1.5\n",
    "y=X*true_w+true_b+np.random.randn(n_samples)*1.0\n",
    "X_train,X_valid,y_train,y_valid=train_test_split(X,y,test_size=0.2)\n",
    "train_size=X_train.shape[0]\n",
    "outlier_size=int(train_size*0.1)\n",
    "outlier=np.array([[10.0,0.0]])+np.random.randn(outlier_size,2)\n",
    "\n",
    "train_data=np.vstack((outlier,np.vstack((X_train,y_train)).transpose()))\n",
    "np.random.shuffle(train_data)\n",
    "\n",
    "X_train=train_data[:,0]\n",
    "y_train=train_data[:,1]\n",
    "\n",
    "plt.scatter(X_train,y_train,marker='o',s=20)\n",
    "X_train=X_train.reshape((-1,1))\n",
    "X_valid=X_valid.reshape((-1,1))\n",
    "\n",
    "def plot_line(w,b,X,c):\n",
    "    y=w*X+b\n",
    "    plt.plot(X,y,linewidth=3.0,color=c)\n",
    "    \n",
    "l1_robust=robust_regression(loss_func=\"l1\")\n",
    "l1_robust.fit(X_train,y_train,X_valid,y_valid)\n",
    "coef=l1_robust.coef()\n",
    "plot_line(coef[\"weight\"],coef[\"bias\"],np.linspace(0.0,10.0,100),c='r')\n",
    "\n",
    "l2_model=robust_regression(loss_func=\"l2\")\n",
    "l2_model.fit(X_train,y_train,X_valid,y_valid)\n",
    "coef=l2_model.coef()\n",
    "plot_line(coef[\"weight\"],coef[\"bias\"],np.linspace(0.0,10.0,100),c='g')\n",
    "\n",
    "l2_model=robust_regression()\n",
    "l2_model.fit(X_train,y_train,X_valid,y_valid)\n",
    "coef=l2_model.coef()\n",
    "plot_line(coef[\"weight\"],coef[\"bias\"],np.linspace(0.0,10.0,100),c='b')\n",
    "\n",
    "plot_line(true_w,true_b,np.linspace(0.0,10.0,100),c='y')\n",
    "plt.legend([\"l1\",\"l2\",\"huber\",\"True\"],fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.Ridge regression\n",
    "One problem with ML estimation is that it can result in overfitting. We can use the MAP estimation to ameliorate this problem.The likelihood function is \n",
    "\\begin{align*}\n",
    "p(D|\\vec{w},\\sigma^2) &=\\prod_{i=1}^N \\mathcal{N}(y_i|\\vec{w}^T\\vec{x_i},\\sigma^2) \\\\\n",
    "                      &=\\left(\\frac{1}{2\\pi\\sigma^2}\\right)^{\\frac{N}{2}}exp \\lbrace -\\frac{\\sum_{i=1}^N (\\vec{w}^T\\vec{x_i}-y_i)^2}{2\\sigma^2} \\rbrace\n",
    "\\end{align*}\n",
    "We can encourage the $\\vec{w}$ to be small, thus resulting in a smoother curve, by using a zero-mean Gaussian prior\n",
    "$$\n",
    "p(\\vec{w})=\\prod_{j=1}^d \\mathcal{N}(w_j|0,\\tau^2)\n",
    "$$\n",
    "Th posterior distribution becomes\n",
    "\\begin{align*}\n",
    "p(\\vec{w}|D) &\\propto p(D|\\vec{w},\\sigma^2)p(\\vec{w}) \\\\\n",
    "             &=\\left(\\frac{1}{2\\pi\\sigma^2}\\right)^{\\frac{N}{2}}exp \\lbrace -\\frac{\\sum_{i=1}^N (\\vec{w}^T\\vec{x_i}-y_i)^2}{2\\sigma^2} \\rbrace    \\prod_{j=1}^d \\mathcal{N}(w_j|0,\\tau^2)\n",
    "\\end{align*}\n",
    "So the **MAP** estimation becomes\n",
    "$$\n",
    "\\hat{\\vec{w}}=\\underset{\\vec{w}}{argmax} \\,p(\\vec{w}|D)\n",
    "$$\n",
    "It's very simple to show that the **MAP** estimation is equivalent to minimizing the following l2 regularization loss\n",
    "$$\n",
    "L(\\vec{w})=\\frac{1}{N}\\sum_{i=1}^N(\\vec{w}^T\\vec{x_i}-y_i)^2+\\lambda \\parallel \\vec{w} \\parallel_2^2\n",
    "$$\n",
    "where $\\lambda$ is a regularization hyperparameter which is chosen by cross validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Bayesian linear regression\n",
    "Although ridge regression is a useful way to compute a point estimation of $\\vec{w}$, but sometimes we want to compute the full posterior over $\\vec{w}$ and $\\sigma^2$\n",
    "### Known $\\sigma^2$\n",
    "The likelihood is given by\n",
    "\\begin{align}\n",
    "p(D|\\vec{w}) &=\\prod_{i=1}^N \\mathcal{N}(y_i|\\vec{w}^T\\vec{x_i},\\sigma^2) \\\\\n",
    "                      & \\propto exp \\left(-\\frac{1}{2\\sigma^2} \\sum_{i=1}^N (y_i-\\vec{w}^T\\vec{x_i})^2 \\right)\n",
    "\\end{align}\n",
    "The conjugate prior to the above Gaussian likelihood is also a Gaussian $p(\\vec{w})=\\mathcal{N}(\\vec{w}|\\vec{w}_0,V_0)$.So the posterior distribution of $\\vec{w}$ is given by\n",
    "\n",
    "$$\n",
    "p(\\vec{w}|D) \\propto \\mathcal{N}(\\vec{w}|\\vec{w}_0,V_0) \\mathcal{N}(D|\\vec{w})=\\mathcal{N}(\\vec{w}|\\vec{w}_N,V_N)\n",
    "$$\n",
    "If $\\vec{w}=\\mathbf{0}$ and $V_0=\\tau^2\\mathbf{I}$, the posterior mean reduces to the ridge estimate.After get the posterior distribution of $\\vec{w}$, we can make prediction according to the distribution.\n",
    "\\begin{align}\n",
    "p(y|\\vec{x},D) &= \\int \\mathcal{N}(y|\\vec{w}^T \\vec{x},\\sigma^2)\\mathcal{N}(\\vec{w}|\\vec{w}_N,V_N)d\\vec{w} \\\\\n",
    "               &=\\mathcal{N}(y|\\vec{w}_N^T\\vec{x},\\sigma_N^2(\\vec{x}))\n",
    "\\end{align}\n",
    "By contrast, the plugin approximation using the **MAP** of the $\\vec{w}$\n",
    "\\begin{align}\n",
    "p(y|\\vec{x},D) &= \\int \\mathcal{N}(y|\\vec{w}^T \\vec{x},\\sigma^2) \\delta_{\\hat{\\vec{w}}} (\\vec{w})d\\vec{w} \\\\\n",
    "               &=\\mathcal{N}(y|\\vec{w}_N^T\\vec{x},\\sigma^2)\n",
    "\\end{align}\n",
    "which under-estimate the uncertainty.\n",
    "### Unknown $\\sigma^2$\n",
    "As usual, the likelihood function has the form\n",
    "\\begin{align}\n",
    "p(D|\\vec{w},\\sigma^2) &=\\prod_{i=1}^N \\mathcal{N}(y_i|\\vec{w}^T\\vec{x_i},\\sigma^2) \\\\\n",
    "                      &=(\\frac{1}{\\sigma^2})^{N/2} exp \\left( -\\frac{1}{2\\sigma^2} \\sum_{i=1}^N (y_i-\\vec{w}^T\\vec{x_i})^2 \\right)\n",
    "\\end{align}\n",
    "We can show that the natural conjugate prior has the following form\n",
    "\\begin{align}\n",
    "p(\\vec{w},\\sigma^2) &=NIG(\\vec{w},\\sigma^2|\\vec{w}_0,V_0,a_0,b_0)  \\\\\n",
    "                    &=\\mathcal{N}(\\vec{w}|\\vec{w}_0,\\sigma^2V_0)IG(\\sigma^2|a_0,b_0) \\\\\n",
    "                    &=\\frac{b_0^{a_0}}{(2\\pi)^{d/2}|V_0|^{\\frac{1}{2}} \\Gamma(a_0)} (\\sigma^2)^{-(a_0+d/2+1)} exp \\left[-\\frac{(\\vec{w}-\\vec{w}_0)^T V_0^{-1}(\\vec{w}-\\vec{w}_0)+2b_0}{2\\sigma^2} \\right]\n",
    "\\end{align}\n",
    "With this prior and likelihood, one can show that the posterior has the following form\n",
    "$$\n",
    "p(\\vec{w},\\sigma^2|D) =NIG(\\vec{w},\\sigma^2|\\vec{w}_N,V_N,a_N,b_N) \n",
    "$$\n",
    "The posterior marginals are as follows:\n",
    "\\begin{align}\n",
    "p(\\sigma^2|D) &=IG(a_N,b_N) \\\\\n",
    "p(\\vec{w}|D)  &=\\mathcal{T}(\\vec{w}_N,\\frac{b_N}{a_N}V_N,2a_N)\n",
    "\\end{align}\n",
    "We can use the posterior distribution of $\\vec{w}$ adn $\\sigma^2$ to compute the posterior prediction which is a **Student T distribution**\n",
    "\\begin{align}\n",
    "p(y|\\vec{x},D) &=\\int \\int \\mathcal{N}(y|\\vec{w}^T\\vec{x},\\sigma^2) p(\\vec{w},\\sigma^2 |D) d\\vec{w}d\\sigma^2 \\\\\n",
    "               &=\\mathcal{T}(y|\\vec{w}_N^T\\vec{x},\\frac{b_N}{a_N}(1+\\vec{x}^T V_N\\vec{x}),2a_N)\n",
    "\\end{align}\n",
    "It's common to set $a_0=b_0=0$, which corresponding to an uninformative prior for $\\sigma^2$, and set $\\vec{w}_0=\\mathbf{0}$ and $V_0=g(X^TX)^{-1}$ for any positive value $g$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:tensorflow-cpu-1.11]",
   "language": "python",
   "name": "conda-env-tensorflow-cpu-1.11-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
